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ABSTRACT 


A new impact debris propagation code has been written to link CTH simulations 
of space debris shield perforation to the Lagrangian finite element code DYNA3D, for 
space structure wall impact simulations. This software (DC3D) simulates debris cloud 
evolution using a nonlinear elastic-plastic deformable particle dynamics model, and 
renders computationally tractable the supercomputer simulation of oblique impacts on 
Whipple shield protected structures. Comparison of three dimensional, oblique impact 
simulations with experimental data shows good agreement over a range of velocities of 
interest in the design of orbital debris shielding. 

Source code developed during this research is provided on the enclosed floppy 
disk. An abstract based on the work described in this report has been submitted to the 
1994 Hypervelocity Impact Symposium. 
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1. Introduction 


The design of the space station Freedom and similar structures for earth orbit 
must include provisions for the effects of hypervelocity impact, which may arise as a 
result of space debris or micrometeorite encounters. This problem can be expected to 
become increasingly important as longer duration space missions are launched, 
increasing the exposure time of orbiting systems. Such missions provide increased 
probability of impact damage while placing greater reliability demands on vehicles and 
structures. Existing light gas gun facilities used in the study of hypervelocity impact 
effects do not generally allow for tests at velocities above ten kilometers per second, 
suggesting the use of computer simulation methods for orbital debris shield design at 
those velocities. 

The accomplishment of design goals for the space station and similar structures 
depends in part upon the development and verification of computationally tractable 
models capable of describing oblique hypervelocity impact effects at velocities beyond 
existing experimental capabilities. Experience to date has shown that the extreme CPU 
time and memory requirements of standard Eulerian hydrocodes (McGlaun et al., 1990) 
make their use in direct simulation of three dimensional impacts on space debris shields 
impractical, even given supercomputer resources. In addition, current Eulerian 
hydrocodes do not in general rigorously account for material history effects on the failure 
of space structures under impact debris loading. Conventional Lagrangian finite element 
codes (Goudreau and Hallquist, 1982), on the other hand, are not suitable for use in 
simulating the shield perforation portion of the impact problem. Mesh distortion effects 
greatly reduce the size of the time step used in the calculations, and mandate frequent 
rezoning. As a result of the preceding difficulties, the only general simulation technique 
demonstrated to date has involved the systematic linking of Eulerian hydrocodes (of 
shield perforation) to Lagrangian finite element models (of debris cloud evolution and 
debris impact on the protected structure). This approach has been developed and 
implemented by the principal investigator (Fahrenthold, 1993). 


The development of particle-based debris cloud evolution models [e.g. Fritts et 
al. (1985), Trease et al. (1990), Monaghan (1988), and Trease (1988)], offers an 
opportunity to directly simulate complete three-dimensional debris impact problems on 
existing supercomputers. Particle-based methods address the previously discussed 
shortcomings of Eulerian and Lagrangian codes in two ways: (1) they effectively 
eliminate the mesh distortion problems of Lagrangian finite element codes, and (2) they 
greatly improve on the CPU time efficiency of Eulerian hydrocodes while allowing for 
accurate tracking of material history dependent effects such as plastic deformation. As an 
application which has severely taxed the capabilities of conventional Eulerian and 
Lagrangian computer codes, orbital debris shield design is well suited to capitalize on the 
strengths of new particle-based modeling techniques. Hence the research presented here 
has: (1) developed a debris cloud evolution code which links shield perforation and wall 
impact simulations, (2) conducted three dimensional, supercomputer based simulations 
of shield impact, debris cloud evolution, and wall impact problems, and (3) evaluated 
and validated the computer models using data from experiments conducted at NASA 
Johnson Space Center. This computer simulation methodology allows for the modeling 
of impacts at velocities beyond ten kilometers per second which are very difficult to 
duplicate in the laboratory. 

2. Methodology 

The problem of hypervelocity impact on space structures has been an object of 
research since the 1950's, as reflected in the summary of Hypervelocity Impact 
Symposia presented during recent conferences in that series (e.g. Anderson, 1986). 
However both experimental and analytical research work has accelerated in recent years, 
with most recent NASA interest focused on debris effects on the space shuttle and 
planned space station. A significant data base exists for impacts at velocities below ten 
kilometers per second (Tower et al., 1987), including studies aimed specifically at NASA 
applications such as space shuttle windows (Schneider and Stilp, 1987) and debris shield 


design (Yew and Kendrick, 1986). However difficulties with conducting experiments at 
higher velocities, relevant to both the space station and future programs, have limited the 
ability to evaluate new designs in a laboratory setting. 

The preceding facts suggest a combined experimental and analytical approach to 
micrometeorite and debris shield design, using experimental data at velocities below ten 
kilometers pier second to critique and verify computer models, which then provide a basis 
for higher velocity design calculations. The simulation work described here was 
conducted on a Cray Y-MP/864 supercomputer at the University of Texas System Center 
for High Performance Computing. 

The principal investigator has made extensive use of existing supercomputer 
based Eulerian and Lagrangian hydrocodes to evaluate their use in orbital debris shield 
design. In general Eulerian codes are best suited to hypervelocity impact simulations 
where impact pressures are sufficiently large to render material strength effects generally 
unimportant. An Eulerian code is best used to predict initial debris cloud mass and 
distribution as a function of projectile material type and velocity, shield material type and 
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geometry, and other parameters. The protected space structure must be designed to avoid 
spallation or fracture under the debris cloud impact. Since spallation and fracture 
processes are very stress and strain history dependent phenomena (Yew and Taylor, 
1992, and Grady and Kipp, 1987), and since Lagrangian hydrocodes are best suited to 
trace stress and strain history in solid materials, a Lagrangian hydrocode model of debris 
cloud impact on the inner wall is most appropriate. The new modeling approach oudined 
here links both parts of the impact event using a new analysis methodology, significantly 
reducing computer resource requirements for oblique impact calculations. 

3. Particle Dynamics Model 

a. Introduction 

This section describes a new modeling approach combining Lagrangian bond 
graphs (Fahrenthold and Wargo, 1994) with a selected finite element discretization 


scheme to allow for direct simulation of the dynamic evolution of debris particles arising 
from hypervelocity impact. 

b. Kinematics 

The homogeneous deformation field associated with the finite element 
discretization employed here allows the position x and velocity v of any mass particle 
"P" in a particular element to be written in the form (Malvern, 1969) 

x = F(t) ( r - r c ) + c(t) ; v = F(t) ( r - r c ) + c(t) (la,b) 

where c(t) and c(t) are the position and velocity of the element center of mass "C", r and 
r c are the position of P and C in the reference (initial undeformed) configuration of the 

body, and F is the deformation gradient tensor. Note that since r and r c in equations (1) 


are constants, the motion of any particle P in the element is determined by the motion of 
C and by the time dependent components of the second order tensor F, related to the rate 
of deformation (D) and velocity gradient (L) tensors by 

D = (1/2) [ L + L T ] ; L = FF 1 (2a,b) 

where the superscripts "-1" and "T" denote the inverse and the transpose. In the special 
case where L is the skew- symmetric tensor whose axial vector is the angular velocity, 
equation (lb) represents rigid body motion (Casey, 1983, and Fahrenthold and Wargo, 
1991a and b). 

c. Kinetic energy 

The homogeneous deformation kinematics of equations (1) allows the kinetic co- 
energy (T ) of a single element of fixed mass "m" and variable volume "V" to be 
expressed in the form 

T* = (1/2) J pv-v dV (3a) 

v 


= d/2) 


{ Jpc*c dV + JpF(r-r c )*F(r-r c )dV + 2c- 
v v 



(3b) 


= (1/2) [ mc*c + tr(F T F J) ] 

where p is the density, "tr" is the trace operator, J is an inertia tensor, 


(3c) 
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J = J p 0 (r-r c )®(r-r c ) dV 0 (3d) 

v 

and conservation of mass requires 

dm = p dV = p Q dV Q (3e) 

with p Q the density in the reference configuration. Note that the third term in equation 

(3b) is zero by definition of the center of mass. Since J is a constant tensor defined in the 
reference configuration, the element momenta p and H are defined by 

p = dT*/()c = me ; H = aT*/9F = FJ (4a, b) 

Since the kinetic energy (T) of the finite element is 

T = p c + tr(H T F ) - T* (5a) 


it follows that 

T = (l/2)[m" 1 p-p + tr{(Hj' 1 ) T H}] (5b) 

The preceding results demonstrate that kinetic energy storage in a single finite 
element may be modeled using the bond graph multiports shown in Figure 1. The kinetic 
energy function may be represented in the most familiar form by the introduction of a 


fourth order inertia tensor G, defined by 
9F9F 

so that 

T* = (1/2) [ m c*c + F:G:F ] ; T = (1/2) [ of 1 pp + H:G -1 :H 
Now the constitutive relation (4b) takes the form 
H = G F 


(6a) 


(6b,c) 

(6d) 


d. Internal energy 

The preceding discussion of inertia effects must be augmented by an RC network 
description of internal energy storage and energy dissipation in the material. Unlike the 
inertia multiports, the capacitance and resistance multiports required to model a particular 
material must be formulated for each material type. To illustrate this procedure, this 
section considers the large strain deformation of an elastic-viscoplastic material (Haupt, 


1985). This case emphasizes the applicability of the bond graph modeling methodology 
described here to very complex engineering materials (Fahrenthold and Wu, 1988), and 
hence to system dynamics and impact dynamics problems of a very general nature. Note 
that the conventional assumptions of infinitesimal strain and linear stress-strain behavior, 
included in the vast majority of mechanical vibrations models, are not adopted here. 
Since the large strains and the relatively complex material response considered here often 
calls for special purpose finite element code development work, the example material 
selected demonstrates the simplicity of the proposed modeling methodology, as 
compared to conventional displacement-based finite element analysis. 

Before formulating a stored energy function for an elastic-plastic material, it is 
appropriate to first define the system kinematics. A very general description of the large 
strain kinematics of elastic-plastic deformation is provided by a multiplicative 
decomposition of the deformation gradient tensor F in the form (Haupt, 1985) 

F = F e FP ; F = F e FP + F e FP (7a,b) 

where F e describes the elastic deformation of the solid material from the unloaded 
plastically deformed configuration and FP relates the unloaded plastically deformed 
configuration of the material to the original undeformed reference configuration. Hence 
the rate of change of F may be decomposed as in equation (7 b) into a first term 
governing the rate of elastic energy storage and a second term governing the rate of 
plastic energy dissipation. In this case, the presence of both complex kinematics and 
energy dissipation effects highlights the value of bond graphs in nonlinear system 
modeling. 

Since most elastic-plastic modeling work begins with a Helmholtz free energy 
function and then derives the internal energy, that procedure is followed here. The 
elastic-plastic kinematics just described are associated with a Helmholtz free energy 
density function of the form (Dashner, 1986) 
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V = \|/(F e , FP, 0) (7c) 

Here is assumed to take the conventional functional form (Bowen, 1989) 

V = (1/po) { (V2) tr(E e ) 2 + p tr(E e2 ) - (3(0-0 o ) tr(E e ) -[(cp a )/(20 o )](0-0 o ) 2 } (7d) 

where X and p are Lame constants for an isotropic solid, c is the specific heat, 0 O is a 
reference temperature, p a is the density in the unloaded plastically deformed 

configuration, and 

E e = (1/2) [ F eT F e -I ] ; p = k (3U2p) ; Po / P(X = det(FP) (7e,f,g) 

with "k" the thermal expansion coefficient and "det" the determinant operator. For an 
element of mass "m", this corresponds to a (conserved) internal energy function 

U = my + 0S (8a) 

where S is the total entropy, defined by the thermodynamic identity 

s = 3(^0 = (m/pa) { cpa f (0/0o) A] + p ^ej } (8b) 

C70 

Equations (7) and (8) may be combined to yield 

U = (m/p Q ) det(FP) { p tr(E e2 ) + (XU) tr(E e ) 2 + P0 O tr(E c ) } + 

[(cm0 o )/2] { [ 1 + S/(mc) - (P/(c Po )) det(FP) tr(E e ) ] 2 - 1 ) (8c) 

The functional form for the internal energy 

U = U(F e , FP, S) (9a) 

leads to the multiport capacitor shown in Figure 2, where in this case 

K e = |^| FP>S = (m/p a ) F*{ [X+((P 2 0 o )/(cp a ))] tr(E e ) I + 

2pE e -[(p0 o S)/(mc)]I} (9b) 

KP = ^ | Fe, S = (m/Pa) ( [(V2)+((p 2 0 o )/(cp a ))] tr(E e ) 2 + 

p tr(E e2 ) - [(P0 o S)/(mc)] tr(E e ) } FP' T (9c) 

0 = | F e F p = e o ( 1+ S/(mc) - [p/(c Pa )] tr(E e ) ) 


(9d) 



e. Plastic deformation 

The materia] description is completed by defining the plastic constitutive relations. 
For simplicity, the rate form 

K d = ti (m/ Po ) LP ; LP = FPfP' 1 (10a,b) 

is adopted here, with ri a viscosity coefficient. The second order tensor K d 
(dimensionally an extensive chemical potential) is the effort power conjugate to LP. 

A model for the elastic-viscoplastic material just defined is shown in Figure 2. 
Note that a transformer with fourth order tensor modulus M, defined by components 


( 11 ) 


^ijrs = d ir^sj 

is introduced in order to conform to the fundamental kinematic relation (10b). 

f. Bond graph model 

The bond graph structure of Figure 2 may be directly augmented with inertia 
multiports representing the kinetic energy contributions associated with p and H. Finally 
if the appropriate RC network model representing elastic-viscoplastic materials is 
introduced, the result is the complete element-level bond graph shown in Figure 3, 
representing the "ith" element of the system. Assuming adiabatic deformation, thermal 
energy is stored, as indicated by the thermomechanical coupling shown in Figure 3. Note 
that in Figure 3 transformers with moduli 

- x. p e(i)_1 . xy|P®_s cPtfH 
^ijkl *jl^ik * ^ijkl “ ^ik^ lj 


(12a, b) 


are introduced in accordance with the kinematics of equations (7a and b). 

g. State equation derivation 

The causally augmented bond graph of Figure 3 and the constitutive relations 
previously defined yield (Rosenberg and Kamopp, 1983) state equations for the "ith" 
element of the form: 

i> (l) = 0 (13a) 


hG) = - MPG) t FPG), S®) 


(13b) 



(13c) 


F e ® = MPO G^)' 1 H® - {p$/Olm®)} MP® M e ® _1 
F p(i) = {p^/CrjmW)} M®* 1 K d (i> (13d) 

S(i> = [l/e( i >(F e ( i ), FP^), S^)] {p^/CTimW)) tr{ K d ^^ T K d ^) } (13e) 


where 

K d (l) = M®' t [ M^)- 7 MP( i ) T K e ( i )(F e ^, FPW, S®) 

- K p ( i >(F e ( i \ FPW, S^))] (13f) 


and the functions K e ®, KP®, and 0® are defined by equations (9). Note that 


M 


e(i)-l 


ijkl = ; MgJ-SikFg'*"; (14a, b,c) 

Equations (13) are nonlinear equations in the state variables: p®, H®, F^*), FPW, and 

SO). 


*P(i) 


w(0-l 


-P(i) 


h. Conclusion 

The outlined modeling methodology may be implemented numerically, for a 
variety of internal energy functions and plasticity models, for use in engineering analysis 
and design. The source code listing at Appendix G represents an isothermal, variable 
compressibility implementation using the plasticity theory of Green and Naghdi (1971). 


4. Comparison to Experiments 

a. Introduction 

Three dimensional simulation of hypervelocity impacts on space structures places 
extreme demands on even supercomputer resources. To reduce the computer time and 
memory requirements of oblique impact simulations, a three-dimensional, deformable 
particle dynamics model of the type just described has been coded and linked to an 
Eulerian hydrocode and a Lagrangian structural code, and applied in the simulation of 
oblique hypervelocity impacts on Whipple shield protected structures. Comparison of the 
results to experimental data shows good agreement at a computer time and memory cost 
much less than that associated with conventional hydrocode calculations. 



This section describes evaluation of the preceding modeling approach using data 
from Whipple shield impact experiments conducted at NASA Johnson Space Center, 
including CPU time requirements for the simulation of representative oblique impact 
problems. 

b. Methodology 

The deformable particle dynamics model of debris cloud evolution is referred to 
here by the title DC3D. This numerical model is used in combination with the Eulerian 
hydrocode CTH (McGlaun et al., 1990) and the structural finite element code DYNA3D 
(Goudreau and Hallquist, 1982) as follows. An Eulerian simulation is first employed to 
model impact on the shield. Then post-processing of the velocity, mass density, and void 
space distribution data from the Eulerian simulation is used to establish the initial state of 
the elastic-plastic model of the debris cloud. Numerical integration of the particle 
dynamics model DC3D, using established system dynamics modeling techniques 
(Rosenberg and Kamopp, 1983) provides a thermodynamically rigorous yet 
computationally efficient basis for predictions of debris cloud evolution over the 
relatively large spaces which normally separate shields from the space structure which 
they protect Output from the debris cloud evolution model then provides a basis for the 
simulation of impact on the wall plate, using DYNA3D. Input data for the wall impact 
simulation is generated automatically by DC3D, based on the final state of the debris 
cloud propagation model. By using a momentum deposition representation of the debris 
cloud loading on the wall plate (calculated from the DC3D results), meshing of each 
debris particle in DYNA3D can be avoided, yielding additional reductions in computer 
time and memory requirements. 

c. Oblique impact example 

The analysis procedure just discussed may be illustrated by considering an 
oblique impact simulation for a Whipple shield configuration. Specifically, consider the 
oblique (23 degree) impact of a 15/64 inch diameter aluminum (2017-T4) sphere, at 7. 1 
kilometers per second, on a 0.063 inch thick aluminum (7075-T6) shield, protecting a 



0.125 inch thick aluminum (7075-T6) wall plate. The shield to wall plate spacing is 4.0 
inches. Figure 4 sh^ws the CTH simulation results for impact on the shield, at four 
microseconds after impact. (Note that due to symmetry, all of the Figures discussed here 
depict only one-half of the physical system modeled.) The code DC3D was then 
employed to: (1) post-processes the CTH results to formulate a debris cloud model, (2) 
propagate the debris to the wall plate, and (3) generate a DYNA3D input file for use in 
modeling impact of the debris cloud on the wall plate. DC3D uses a variable step method 
for numerical integration of the first order state equations describing the elastic-plastic 
debris cloud dynamics, and commercial plotting routines for graphical representation of 
the simulation results. Figure 5 shows the predicted impact momentum distribution on 
the wall plate. Finally a DYNA3D simulation of the wall plate impact was conducted, 
using a momentum deposition approach to represent the effects of the debris. Figures 6 
and 7 show front and rear views of the wall plate impact simulation, including extensive 
spallation observed in the experiment 

By making appropriate use of Eulerian and Lagrangian codes for those pans of 
the simulation where they are both accurate and computationally tractable, this approach 
does not suffer from the limiting assumptions of many previous attempts to model shield 
impact problems (e.g. Swift et al., 1983). In addition, that portion of the total impact 
simulation described by the debris cloud model incorporates arbitrarily nonuniform, 
three-dimensional velocity, density, and void space distributions not admitted by many 
simplified debris models published to date (Grady and Passman, 1990). In summary, the 
approach used here makes use of the known strengths of available codes while providing 
an essential improvement in computational efficiency for oblique impact simulation. Such 
an improvement is essential before computer codes can provide a practical analytical 
design tool. 

The preceding example demonstrates the feasibility of the modeling approach. Its 
computational efficiency is such that computer resource requirements are relatively 



modest, considering the three dimensional nature of the simulations. Typical CPU time 
requirements are listed in Table 1. 

d. Ballistic limit calculations 

Finally it is worthwhile to compare some hypervelocity impact simulation results 
with ballistic limit calculations based on published experimental research. Figure 8 
depicts a typical Ballistic limit curve, and indicates the velocities and impact panicle 
diameters used in the simulations. Figures 9 through 20 show results of the simulation of 
oblique Whipple shield impacts at 7 and 10 kilometers per second, using the debris 
propagation code and the general modeling methodology developed under this project. 
The cases represented are listed in Table 2. 

Note that the simulations were conducted for particle diameters above and below 
the expected ballistic limits for the two velocities. The simulations are consistent with the 
ballistic limit curve data at 7 km/sec. At 10 km/sec, the simulations suggest slightly less 
wall damage than might be expected from the ballistic limit curve. Computer time 
requirements are very modest, considering the complexity of the problems. 

5. Conclusion 

The outlined research has addressed fundamental problems relevant to the 
development of space vehicles and structures for a variety of missions. It makes use of 
state of the art supercomputer resources while applying the newest numerical modeling 
techniques, designed to make oblique impact simulations computationally tractable. It 
included numerical implementation of a nonlinear, elastic -plastic debris cloud dynamics 
model important to accurate shield impact calculations. The resulting simulation capability 
can provide an important adjunct to experimental work on space shield design for a 
variety of future missions. The outlined work has been coordinated directly with NASA 
JSC experimental research efforts, and the final source code has been provided to 
NASA. It is therefore suggested that the proposed research has addressed important 
objectives of the NASA Regional Universities Grant Progr am 
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Table 1. Example Oblique Impact Simulations 


1 5 


shield to wall spacing: 
shield thickness: 
wall thickness: 
impact velocity: 
impact obliquity: 
impact mass: 
material: 


4.0- 8.0 inches 
0.063-0.071 inches 
0.125-0.160 inches 

7. 1- 7.6 km/sec 
23.0-58.4 degrees 
0.202-0.376 grams 
aluminum 


Simulation 

type 

Computer 

system 

Code 

Memory 

CPU time 

shield impact 

Cray Y-MP/864 

CTH 

<64 MB 

2.0-6.0 hrs 

debris evolution 

IBM RS/6000* 

DC3D 

< 16 MB 

0.5-3.0 days 

wall impact 

Cray Y-MP/864 

DYNA3D 

<64 MB 

0.5- 1.0 hrs 


*The CPU time range given is for a Model 320, a very low performance system. 



Table 2. Ballistic 



bumper thickness 
wall thickness 
impact obliquity 
material 


(also see the ballistic 

Simulation # 

Impact velocity 

7a 

7 km/sec 

7b 

7 km/sec 

10a 

10 km/sec 

10b 

10 km/sec 


bumper impact simulations (CTH): 
debris transport calculations (DC3D): 
wall impact calculations (DYNA3D): 


imit Simulations 
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= 0.127 cm 

= 0.3175 cm 

= 25 degrees 

= A1 6061-T6 

imit plot in Figure 8) 


0.575 cm 
0.475 cm 
0.530 cm 
0.430 cm 


Figures 

9 through 1 1 
12 through 14 
15 through 17 
18 through 20 


1.80-2.11 CPU hrs, Cray YMP 
3.34-5.67 CPU hrs, IBM RS/6000 
0.27-0.28 CPU hrs, Cray YMP 


List of Figures 


Figure 1. Inertia multiports: deformable particle dynamics model. 

Figure 2. RC multiports: deformable particle dynamics model. 

Figure 3. Bond graph: deformable particle dynamics model. 

Figure 4. CTH simulation of oblique impact on a Whipple shield (t = 4 (isec). 

Figure 5. DC3D simulation results for impact momentum distribution on the wall plate. 

Figure 6. DYNA3D simulation results: front surface of the wall plate. 

Figure 7. DYNA3D simulation results: rear surface of the wall plate. 

Figure 8. Ballistic limit plot. 

Figure 9. CTH simulation 7a results: the shield and the debris cloud at four 
microseconds after impact. 

Figure 10. DC3D simulation 7a results: areal density of the debris cloud momentum 
deposited on the wall plate. 

Figure 11. DYNA3D simulation 7a results: front surface of the wall plate. 

Figure 12. CTH simulation 7b results: the shield and the debris cloud at four 
microseconds after impact. 

Figure 13. DC3D simulation 7b results: areal density of the debris cloud momentum 
deposited on the wall plate. 

Figure 14. DYNA3D simulation 7b results: front surface of the wall plate. 

Figtire 15. CTH simulation 10a results: the shield and the debris cloud at three 
microseconds after impact. 

Figure 16. DC3D simulation 10a results: areal density of the debris cloud momentum 
deposited on the wall plate. 

Figure 17. DYNA3D simulation 10a results: front surface of the wall plate. 

Figure 18. CTH simulation 10b results: the shield and the debris cloud at three 
microseconds after impact. 

Figure 19. DC3D simulation 10b results: areal density of the debris cloud momentum 
deposited on the wall plate. 

Figure 20. DYNA3D simulation 10b results: front surface of the wall plate. 
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APPENDIX A 



Example Input File (CTH simulation) 


] 
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* 

*eor* cgenin 
* 

* 

* cthgen input 

* 

* 

* Title record 

* 

Oblique whipple shield impact (7a) 
* 

* 

* control records 

* 

control 

ep 

* mmp 

endc 

* 

* 

* edit records 

* 

edit 

block 1 
expanded 
endb 
ende 
* 


* mesh records 

* 

mesh 

block 1 geom=3dr type~e 


x0 

0.0 



xl 

n=50 

w— 1.5875 

rat-1 

endx 




yO - 

11.5 



yi 

n=125 

w=3. 96875 

rat=l 

endy 




zO 

- 1.000 



zl 

n=125 

w-3 .96875 

rat=l 


endz 

xact - 0.0 1.00 

yact = -11.5 -10.40 

zact - -0.50 1.00 

endb 
endm 
* 

*********************************************************************** 

* 

* material insertion records 
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insertion of material 
* 

block 1 

* 

package projectile 
material 1 
numsub 49 

velocities xvel 0.0 yvel 6.3441545e5 zvel 2.9583278e5 
insert sphere 

center 0.0 -10.5745 0.0 

r 0.2875 

endi 
endp 

* 

package shield 
material 1 
numsub 49 

velocities xvel 0. yvel 0. zvel 0. 
insert box 

xl 0.0 yl -10.287 zl -1.00 x2 1.5875 y2 -10.16 z2 2.96875 
endi 
endp 

* 

endb 

* 

endi 

* 

*********************************************************************** 

* 

* eos records 

* 

* eos num 2 
eos 

matl mgrun eos=6061-t 6_al 

* matl mgrun eos=2024-t 4_al 

* mat2 mgrun eos=7075-t 6_al 

* 

* aneosl -1 'Aluminum library* lib=6 type=4 rhug=-l. thug=-l . 

* 

ende 

* 

*********************************************************************** 

* 

* material strength records 

* 


epdata 
* matep 

1 

jf rac=7 

* matep 

2 

jf rac=8 

matep 

1 

st=6061-t6 aluminum 

* matep 

1 

st=202 4-t 4_aluminum 

* matep 

2 

st=7075-t 6_aluminum 

* matep 

1 

st=user r0=1.0 

* matep 

2 

st=user rO^l.O 

mix = 

5 


ende 

* 




*********************************************************************** 

* 
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endinput 

* 

* end of cthgen input 

* 

* 

*eor* cthin 
* 

* 

* cth input 

* 

* 

* Title record 

* 

Oblique whipple shield impact (7a) 
* 

* 

* control records 

* 


restart 
nu “ 1 
endr 
* 

control 

t st - 1.0e-6 
nsc = 5000 
cpshift =60. 

* mmp 
endc 

* 

convect 

convection = 1 
endc 
★ 

*********************************************************************** 

* 

*_ time step records 
* 

mindt 

time =0. dt = l.e-12 
endn 
* 

maxdt 

t ime=0 . dt=0 . 5e-10 

time =3 . e-10 dt=l . 
endx 
* 


* 


* tracer records 

* 


tracer 

* add 0 . 0 

* add 0 . 0 

* add 0 . 0 


-10.528 0.0 
-10.368 0.0 
-10.287 0.0 
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* 

add 

0.0 

-9.368 





* 

add 

0.0 

-0 . 127 





* 

add 

0.0 

-9.368 

to 

1,0 

-9.368 

n 3 

* 

add 

0.0 

-0 . 127 

to 

1.0 

-0 . 127 

n 9 


endt 


* 

* 

* edit records 

* 




* * 


edit 
shortt 
t ime=0 . 
t ime=l . e-5 
time=l .e-4 


dtf =5 . e-6 
dtf =1 . e-5 
dtf=l . e-4 


ends 
longt 
time=0 . 
time=l . e-5 
time=l .e-4 
endl 
plott 


dtf =5 . e-6 
dtf =1 . e-5 
dtf =1 .e-4 


time=0 . dtf=0.5e-6 
time=15.0e-6 -dtf=l . Oe-6 
time=30.0e-6 dtf=2.0e-6 


endp 

histt 



t ime = 0 , 

dtf = 1. 

e-8 

* 

ht racerl 




* 

htracer2 




* 

htracer3 




* 

bxyz 1 

0.0 

-8.01 

0.0 

★ 

bxyz 1 

0.2 

-8.01 

0.0 

* 

bxyz 1 

0.4 

-8.01 

0.0 

* 

bxyz 1 

0.6 

-8.01 

0.0 

* 

bxyz 1 

0.8 

-8.01 

0.0 

★ 

bxyz 1 

1.0 

-8.01 

0.0 


endh 

ende 


* 


* boundary condition records 

* 

boundary 
bhydro 
block 1 

bxb=0 bxt=l byb=l byt=l bzb=l bzt=l 
endb 
endh 
endb 
* 

* end of cth input 

* 

*********************************************************************** 


* 

*eor* hisinp 
* 
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* hisplt input 

* 

*********************************************************************** 

* 

bottom=of f 

plot time cpu v2=dt v3=etot 

*plot time pressure. 1 v2=pressure . 2 v3=pressure . 3 
*plot time yvelocity.l v2=yvelocity , 2 v3=yvelocity . 3 
*plot time zvelocity.l v2=zvelocity . 2 v3=zvelocity . 3 
*plot time pressure. 4 v2=pressure . 5 v3=pressure . 6 
*plot time pressure. 7 v2=pressure . 8 v3=pressure . 9 


* end of hisplt input 

* 

* 

*eor* pltinp 
* 

* 

* cthplt input 

* 


* 

units cgsk 
bottom=of f 
* 

*color table - 6 

*color material - 16 112 

*color nmaterial * 31 127 

*3dplot 

*tcolor = 6 

* 

*3dline, mat=0 
3dortho, mat=0 
* 

*2df ix x=0 . 0 

*2dplot if tracer dots=density=3 . 0 
- * 

* end of cthplt input 

* 
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~ APPENDIX B 


Example Input File (CTH simulation restart) 


r i 


hd 
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A********************************************************************** 

* 

*eor* cgenin 
* 

*********************************************************************** 

* 

* cthgen input 

* 

*********************************************************************** 

* 

* Title record 

* 

Oblique whipple shield impact (7a) 

* 

*********************************************************************** 

* 

* control records 

* 

control 

ep 

* mmp 

endc 

* 

*********************************************************************** 

* 

* edit records 

* 

edit 

block 1 
expanded 
endb 
ende 
* 

*********************************************************************** 

* 

* mesh records 

* 

mesh 

block 1 geom=3dr type=e 
xO 0.0 

xl n=50 w*1.5875 rat-1 . . 

endx 

yO -11.5 

yl n=125 w=3 . 96875 rat=l . 

endy 

zO -1.000 


zl n ; 

-125 w=3. 

96875 rat-1. 

endz 
xact * 

0.0 

1.00 

yact = 

-11.5 

-10.40 

zact * 

-0.50 

1.00 

endb 



endm 

+ 



*********************************************************************** 

* material 

* 

insertion 

records 
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insertion of material 
* 

block 1 

* 

package projectile 
material 1 
numsub 49 

velocities xvel 0.0 yvel 6.3441545e5 zvel 2.9583278e5 
insert sphere 

center 0.0 -10.5745 0.0 

r 0.2875 

endi 
endp 

* 

package shield 
material 1 
numsub 49 

velocities xvel 0. yvel 0. zvel 0. 
insert box 

xl 0.0 yl -10.287 zl -1.00 x2 1.5875 y2 -10.16 z2 2.96875 
endi 
endp 

* 

endb 

★ 

*endi 

* 

*************************************************************** ******** 
* 

* eos records 

* 

* eos num 2 
eos 

mat 1 mgrun eos=6061-t 6_al 

* matl mgrun eos-2024-t 4_al 

* mat2 mgrun eos=7075-t6_al 

* 

* aneosl -1 ’Aluminum library 1 lib=6 type =4 rhug=-l. thug=-l . 

* 

ende 

* 

*********************************************************************** 

★ 

* material strength records 

* 

epdata 

* matep 1 jfrac=7 

* matep 2 jfrac=8 

matep 1 st=6061-t 6_aluminum 

* matep 1 st=2024-t 4_aluminum 

* matep 2 st=7075-t 6_aluminum 

* matep 1 st=user r0=1.0 

* matep 2 st=user r0=1.0 
mix * 5 

ende 

* 

*********************************************************************** 


* 
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endinput 

* 

* end of cthgen input 

* 

*********************************************************************** 

* 

*eor* cthin 

* 

*********************************************************************** 

* 

* cth input 

* 

*********************************************************************** 

* 

* Title record 

* 

Oblique whipple shield impact (7a) 

* 

*********************************************************************** 

* 

* control records 

* 

restart 
newfile 
nu = 3 
endr 
* 

control 

tst = 4.0e-6 
nsc = 5000 
cpshift = 60 . 

* mmp 
endc 

* 

convect 

convection =* 1 

endc 

* 

*********************************************************************** 

* 

* time step records 

* 

mindt 

time =0. dt = l.e-12 
endn 
* 

maxdt 

time=0 . dt=0 . 5e-10 

time =3 . e-10 dt=l . 
endx 
* 

*********************************************************************** 

* 

* tracer records 

* 


tracer 

* add 0.0 -10.528 0.0 

* add 0.0 -10.368 0.0 
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* 

add 

0.0 

-10.287 

o 

o 




* 

add 

0.0 

-9.368 





* 

add 

0.0 

-0 . 127 





* 

add 

0 . 0 

-9.368 

to 

1.0 

-9.368 

n 3 

* 

add 

0.0 

-0 . 127 

to 

1.0 

-0 . 127 

n 9 


endt 

* 

* 

* edit records 

* 


edit 
shortt 
t ime=0 . 
time=l . e-5 
t ime=l . e-4 
ends 
longt 
time=0 . 
time=l . e-5 
t ime=l . e-4 
endl 
plott 
t ime=0 . 
time=15 . Oe- 
time=30 . Oe- 
endp 
histt 


dt f =5 . e-6 
dt f =1 .e-5 
dtf=l .e-4 


dtf=5 .e-6 
dtf=l .e-5 
dtf=l .e-4 


dtf =1 . 5e-6 
•6 dtf =1 . Oe-6 
-6 dtf =2 . Oe-6 



t ime = 0 . dt f » 1 , 

. e-8 

* 

ht racerl 


* 

ht racer2 


* 

htracer3 


* 

bxyz 1 0.0 -8.01 

0.0 

* 

bxyz 1 0.2 -8.01 

0.0 

* 

bxyz 1 0.4 -8.01 

0.0 

* 

bxyz 1 0.6-8.01 

0.0 

* 

bxyz 1 0.8 -8.01 

0.0 

* 

bxyz 1 1.0 -8.01 

0.0 


endh 

_ende 

* 


* boundary condition records 

* 

boundary 
bhydro 
block 1 

bxb-0 bxt=l byb=l byt=l bzb=l bzt=l 
endb 
endh 
endb 
* 

* end of cth input 

* 

***************************************************************** ****** 


*eor* hisinp 
* 
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* 

* hisplt input 

* 

* 


bottom=of f 
plot time cpu 


*plot time 
*plot time 
*plot time 
*plot time 
*plot time 


pressure . 1 
yvelocity . 
zvelocity . 
pressure . 4 
pressure . 7 


v2=dt 

v2=pressure . 2 
1 v2=yvelocity . 2 
1 v2=zvelocity . 2 
v2=pressure . 5 
v2 =s pressure . 8 


v3=etot 
v3= s pressure . 3 
v3=yvelocity . 
v3=zvelocity . 
v3=pressure . 6 
vS^pressure . 9 


3 

3 


* 


* end of hisplt input 

* 




* 


*eor* pltinp 
* 


* 

* cthplt input 

* 


*********************************************************************** 

* 

units cgsk 
bottom=of f 
* 

*color table = 6 

*color material * 16 112 

*color nmaterial * 31 127 

*3dplot 

*tcolor = 6 

* 

*3dline, mat~0 

3dortho, mat=0 
* 

*2d£ix x=0 . 0 

*2dplot if tracer dots=density=3 . 0 
* 

* end of cthplt input 

* 

************************************************************************ 
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Example Input File (DC3D) 
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* dc3d input file 

* 


* 

itype 

jtype 






2 

1 





* 

il 

i2 






2 

33 





★ 

jl 

j2 






64 

112 





* 

kl 

k2 






17 

80 





* 

xlow 


delta x 





0.0 


0.031750 





ylow 


delta y 





-11.5000 


0.031750 




* 

zlow 


delta z 





- 1.0000 


0.031750 




* 

tout 


tmin 

dencut 




25.0 


0.0 

0.01 



* 

idmax 

jdmax 

kdmax 

idmaxn jdmaxn 

kdmaxn 


51 

5 

101 

51 

5 

101 

* 

xdmin 


xdmax 

xdminn 


xdmaxn 


0.0 


8.00 

0.0 


4.00 

* 

ydmin 


ydmax 

ydminn 


ydmaxn 


0.0 


0.3175 

0.0 


0.3175 

* 

zdmin 


zdmax 

zdminn 


zdmaxn 


-4.00 


12.00 

2.0 


10.0 

* 

vmu 


vlamb 

denref 


eta 


0.259 


0.503 

2.703 


1.0e-3 

* 

zeta 


eps 

ymxref 




1.0e-9 


1.0e-6 

1 . Oe-6 



* 

mdf actor 







1.0 








D-l 


APPENDIX D 

Example Input File Header (DYNA3D) 
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Debris Cloud 

Wall Impact Model 




88 

2 25763 20001 





0 0 

0 0 0 0 0 

0 

0 0 




e20 . 0 





0 


1 

0 

0 

0 1619 

100.0e-01 

-2 . e+4 50 . e-01 


1.0e-18 



0 

0 

0 

1 13 

2.703 






e-p w/failure 

0.276 2 . 900e-03 3.88e-03 1.500e-0 -12.et33 

7 . 20e-l 


2 1 1.0 
elastic 

1.0 

0.3 


APPENDIX E 


Plotting Routine Source Code (DCMPLT) 
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c program dcmplt 

c 

dimension x (500) , y (500) , z (500) , 

1 zz(500, 500) , cval (10) 
c 

external grctr,efsplt 
c external grct r , egsgl, ef spit 

c 

open ( 1 , file- 1 depost * pit 1 ) 
c open (2, f ile= T dc3din * ) 

c open ( 9 , f i le= 1 dcmplt . out ' ) 

c 

c do 11 i=l,500 

c x ( i) =0 . 0 

c y (i) =0 . 0 

c z ( i) =0 . 0 

c 11 continue 
c 

c scale=l . Oe+6 

scale=l . Oe+3 
c 

jmax=l 
kmax=l 
zmax=0 . 0 
c 

do 10 i-1, 500000 

read ( 1, 101, end=98) j , k, xx, yy, pmag 
101 format (2i6, 3el5 . 3) 
x ( j ) =xx 
y (k) =yy 

zz ( j , k) =pmag*scale 
if (zz(j,k) .gt.zmax) zmax=zz(j f k) 
if ( j .gt . jmax) jmax= j 
if (k . gt . kmax) kmax-k 
10 continue 
c 

9 8 dummy = 1 . 0 
c 

iunit=0 

ldzz=500 

c iopt-1+2+4 

iopt=2+8 
ncv=8 
c 

cval ( 1) **1 . 0e-33 
cval (2) =zmax 
c 

call gretr ( jmax, kmax, x, y, zz, ldzz, iopt , nev, cval) 
call egsgl ( 1 ! contour_j?lot use$ ' , 1 dcmplt .dl$ 1 ) 

call egsgl ( 1 . 1 use$ T , ' dcmplt .d2$ 1 ) 

call egsgl ( 1 . 1 viewport$ f , 0 . 1, 0 . 9, 0 . 1, 0 . 9) 
call ef spit (iunit , 1 1 ) 


c 


stop 

end 
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Post-processor Source Code (DCPOST) 
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L 




program depost 
c 

c**** program depost 
c 

dimension if lg (100000) , px (100000) , py ( 100000 ) , pz (100000 ) , 
1 jenew (100000) 

dimension xcen (500) , zeen (500) , emom(500, 500) 
c 

character* 8 chO, chi, ch2, ch3, ch4, ch5, ch6, ch7, ch8, ch9 
c 

open(l,file ^dcSdin 1 ) 
open ( 12 , f ile^ 1 dc3d . dyn T ) 
open ( 6 1 , f i le= 1 dc3d . mom f ) 
c open (13, file= 'debase 1 ) 

c 

open (21, f ile= 'depost .out ' ) 
open (22 , f ile= ' depost . pit 1 ) 
c 

c time=l . 0 

time=0 . 0 
c 

c factor=0.75 

c factor=1.0 

c 

read(l,110) itype, jtype 

110 format (///2il0) 

read(l, HI) il, i2, jl, j2, kl, k2 

111 format ( /2il0) 

read ( 1, 112 ) xlow, deltx, ylow, deity, zlow, deltz 

112 format (/2el5 . 6) 
read(l, 113) tout 

113 format (/el5 . 3) 

read (1,114) idmax, jdmax, kdmax, idmaxn, jdmaxn, kdmaxn 

114 format (/6il0) 

read (1, 115) xdmin, xdmax, xdminn, xdmaxn, 

1 ydmin, ydmax, ydminn, ydmaxn, 

2 zdmin, zdmax, zdminn, zdmaxn 

115 format (/4el5 . 3) 
read(l, 116) factor 

116 format (/////e 15.3) 
c 

nent 331 idmax* jdmax* kdmaxt (idmax-1) * ( jdmax-1) * (kdmax-1) + 

1 8+1+5 


c if ( jtype . ne . 1) go to 777 

if ( jtype . eq. 0 ) go to 111 
c 

do 14 i=l,ncnt 

read (12, 202) chO, chi, ch2 , ch3, ch4, ch5, ch6, ch7, ch8 , ch9 

202 format (10a8) 
c 

if (jtype. eq. 3) go to 14 
c 

write (21,203) ch0,chl, ch2, ch3, ch4, ch5, ch6, ch7, ch8, ch9 

203 format (10a8) 

14 continue 


111 dummy=l . 0 
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do 10 j = l, 100000 
iflg( j)=0 
px ( j) =0 . 0 
py ( j ) -0 . 0 
pz ( j) =0 . 0 

10 continue 
c 

do 11 i-1, 500000 

c read (12, 102,end=75) j , dpx, dpy, dpz 

read (61, 102,end=75) j , dpx, dpy, dpz 
102 format (i8, 3el0 . 3) 
if lg ( j) =1' 

px ( j ) =px ( j ) +dpx*f actor 
py ( j ) =py ( j ) +dpy*f actor 
pz ( j ) =pz ( j ) +dpz* factor 

11 continue 
c 

75 dummy =1 . 0 
c 

r i dma x= i dma x 
rkdmax^kdmax 

deldx= (xdmax-xdmin) / (ridmax-1 . 0) 
deldz= ( zdmax-zdmin) / ( rkdmax-1 . 0 ) 
r i dmxn= idma xn 
rkdmxn=kdmaxn 

deldxn= (xdmaxn- xdminn) / (ridmxn-1 . 0) 
deldzn= ( zdmaxn-zdminn) / (rkdmxn-1 . 0) 
c 

do 710 k=l, kdmax-1 
do 711 i=l,idmax-l 
rk=k 
ri=i 

zcen (k) *zdmin+deldz/2 . 0+ ( rk-1 . 0 ) *deldz 
xcen (i) =xdmin+deldx/2 . 0+ (ri-1 . 0) *deldx 
emom(i, k) =0 . 0 
c 

jelem=i+ (idmax-1) * (k-1) 
jenew ( jelem) =0 

if (xcen (i) . It . xdminn) go to 711 
if (xcen (i) .gt .xdmaxn) go to 711 
if (zcen (k) . It . zdminn) go to 711 
if (zcen (k) . gt • zdmaxn) go to 711 
inew-int ( (xcen (i) -xdminn) /deldxn) +1 
knew=int ( (zcen (k) -zdminn) /deldzn) +1 
jenew ( jelem) =inew+ (idmaxn-1) * (knew-1) 
c 

711 continue 
710 continue 
c 

do 12 j=l, 100000 
c 

j out — j 

if ( jtype.eq.3) jout= jenew ( j ) 
if ( jout .eq. 0) go to 12 
c 

if (if lg ( j) .eq. 0) go to 12 

write (21,201) jout, px ( j ) , py ( j) ,pz ( j) r time 
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201 format (18, 4e 10.3) 
c 

if ( jtype.eq.3) go to 12 
c 

iel=mod ( j , idmax-1) 
kel=j/ (idmax-1) +1 

pmag= (px ( j ) **2+py (j) **2+pz { j ) **2 ) **0,5 
emom (iel, kel) =pmag/ (deldx*deldz) 
c emom ( iel, kel) =pmag 

c riel^iel 

c rkel=kel 

c xi=xdmin+deldx+ (riel-1 . 0 ) *deldx 

c zi=zdmin+deldz+ (rkel-1 , 0 ) *deldz 

c write (22 , 212 ) iel, kel, xi, zi, pmag, time 

c 212 format (2i6, 4el5 . 3) 
c 

12 continue 
c 

if ( jtype . eq. 3) go to 99 
c 

do 721 k=l,kdmax-l 
do 722 i»l, idmax-1 

write (22,213) i,k, xcen (i) , zcen (k) , emom(i, k) 
213 f ormat (2i6, 3el5 . 5) 

722 continue 
721 continue 
c 

9 9 dummy =2 . 0 
c 

stop 

end 
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program dc3d 


c**** program dc3d 



common/atype/itype, jtype 

common/mesh/il, i2, j 1, j2 , kl , k2 , deltx, deity, deltz , 

■ - : 


1 xlow, ylow, zlow 

r- ; 


common / wa Ida t/ywa 11, tcut , tmin, dencut , 



1 idmax, jdmax, kdmax, xdmin, xdmax, ydmin, ydmax, zdmin, zdmax 
common/walnew/idmaxn, jdmaxn, kdmaxn, 

: ! 5 


1 xdminn, xdmaxn, ydminn, ydmaxn, zdminn, zdmaxn 

m ■ 


common /props / den, denref , rmass, vmu, vlamb, eta, zeta, phi, r j (3, 3 ) 



common/ intpar/epsint , ymxref 

. 

c 

open (1, f ile= T dc3din T ) 


c 

open ( 7 , f ile= 1 pltdat / dc3d . pit 1 1 ) 



open ( 8 , f ile= 1 pltdat /dc3d . pit 2 1 ) 

: - : 


open (4, f ile= ’dc3d. cth T ) 

— 


open ( 9 , f ile= ' dc3d . dyn ' ) 


c 

open ( 61 , f ile= ! dc3d .mom’ ) 

— 

c 

open (10 , f ile= 1 dc3d. dbg T ) 


c 

read(l,110) itype, jtype 
110 format (///2il0) 



read(l. Ill) il, i2, jl, j2, kl*, k2 
111 format (/2il0) 

read(l, 112) xlow, deltx, ylow, deity, zlow, deltz 



112 format ( /2el5 . 6) 

read ( 1 , 113) tcut , tmin, dencut 

113 format ( /3el5 . 3) 

read (1, 114) idmax, jdmax, kdmax, idmaxn, jdmaxn, kdmaxn 

— 


114 f ormat ( / 6il0 ) 

read (1,115) xdmin, xdmax, xdminn, xdmaxn. 



1 ydmin, ydmax, ydminn, ydmaxn, 

2 zdmin, zdmax, zdminn, zdmaxn 



115 format </4el5 . 3) 



read(l,116) vmu, vlamb, denref , eta 

- - 


116 format (/4el5 . 3) 



read(l, 117) zeta, epsint, ymxref 



117 f ormat ( /3el5 . 3) 

fu 

c 

ywall=ydmin 


c 



c 

call rdavs 


c 

c 

if (itype .eq. 2) call dyngen 



if (itype . eq. 2 . and. jtype . eq. 1) call dyngen 
if (itype .eq. 2 . and . jtype . eq. 2) call dyngen 
if (itype.eq.2 .and. jtype.eq.3) call dyngen 
if (jtype. eq. 2) go to 99 

w 


if ( jtype .eq. 3) go to 99 


c 

call rdavs 

Ww 

c 

if st=2 
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I 1 


f= 

e 3 

1 J 


u 


jf st=2 
kf st=2 
c 

c ilast=2 

c jlast=2 

c klast=2 

c 

ilast=12-il 
jlast= j2- jl 
klast=k2-kl 
c 

lfst»l 

llast^S 

c 

do 10 i=ifst,ilast 
do 11 jfst , jlast 
do 12 k=kfst,klast 
do 13 l=lfst,llast 
call dcelem (i, j , k, 1) 

13 continue 
12 continue 
11 continue 
10 continue 
c 

9 9 dummy= 1 . 0 
c 

stop 

end 

c 

subroutine dyngen 
c 

c**** subroutine dyngen 
c 

common/ a type/itype, jtype 

common /walnew/ idmaxn, jdmaxn, kdmaxn, 

1 xdminn, xdmaxn, ydminn, ydmaxn, zdminn, zdmaxn 
common/waldat /ywall, tcut , tmin, dencut, 

1 idmax, jdmax, kdmax, xdmin, xdmax, ydmin, ydmax, zdmin, zdmax 
c 

dimension x (8) ,y (8) , z (8) , xn (501) ,yn<501) ,zn(501) 
c 

if (jtype. ne. 3) go to 88 
idmax=idmaxn 
jdmax= jdmaxn 
kdmax=kdmaxn 
xdmin=xdminn 
ydmin=ydminn 
zdmin= zdminn 
xdmax= xdmaxn 
ydmax=ydmaxn 
zdmax= zdmaxn 
88 dummy-1.0 
c 

r idmax- idmax 

delx= (xdmax-xdmin) / (ridmax-1 .0) 
r jdmax= jdmax 

dely= (ydmax-ydmin) / (r jdmax- 1 .0) 
rkdmax=kdmax 
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delz= (zdmax-zdmin) / (rkdmax-1 . 0) 
c 

do 10 i=l, idmax 
ri=i 

xn(i)=(ri-1.0) *delx 

10 continue 
c 

do 11 j=l , jdmax 
r j = j 

yn(j)=(rj-1.0) *dely 

11 continue 
c 

do 12 k=l f kdmax 
rk=k 

zn (k) = (rk-1 .0) *delz 

12 continue 
c 

c define the nodes 

c 

c nconst=0 

c 

do 20 j=l, jdmax 
do 21 k=l, kdmax 
do 22 i=l, idmax 
c 

nconst=0 

if (i . eq. 1) nconst=l 

c 

go to 761 
c 

if (i .eq. idmax) nconst=4 
if (k.eq. 1) nconst-5 

if (k.eq. kdmax) nconst^S 
if (i . eq. 1 . and. k . eq. 1) nconst-7 

if ( i . eq . 1 . and .k.eq. kdmax) nconst=7 

if (i .eq. idmax . and . k . eq . 1) neons t=7 

if ( i . eq . idmax . and .k.eq. kdmax) nconat=7 
c 

761 dummy- 1 . 0 
c 

nnode=i+idmax* (k-1) + idmax* kdmax* ( j-1) 
write (9, 101) nnode f nconst, xn (i) ,yn( j) /zn(k) 
c 101 format (i5, i5 r 3e20 . 8) 

101 format (i8, i5 f 3e20 . 8) 

22 continue 
21 continue 
20 continue 
c 

nconst=*0 

c 

x { 1) =xdmin 
y ( 1 ) =ydmin- 1 . 0 
z ( 1) =zdmin 
x (2) =xdmin+l . 0 
y (2) =ydmin-l . 0 
z (2) =zdmin 
x (3) =xdmin+l . 0 
y (3) =ydmin-l . 0 
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z (3) =zdmin+l . 0 
x (4) =xdmin 
y ( 4) =ydmin-l . 0 
z ( 4 ) =zdmin+l . 0 
x (5) =xdxnin 
y (5) =ydmin-2 . 0 
z { 5) =zdmin 
x ( 6 ) =xdmin+l . 0 
y ( 6) =ydmin-2 . 0 
z ( 6) =zdmin 
x (7) =xdmin+l . 0 
y (7 ) =ydmin-2 . 0 
z (7 ) =zdmin+l . 0 
x ( 8) =xdmin 
y ( 8) =ydmin-2 . 0 
z ( 8 ) =zdmin+l . 0 
c 

do 40 1=1,8 

nnode=idmax+ idmax* (kdmax-1) 4-idmax*kdmax* ( jdmax-1) 4-1 
write (9,103) nnode, neons t , x(l) ,y(l) ,z(l) 
c 103 format (i5, i5, 3e20 . 8) 

103 format (18,15, 3e20 . 8) 

40 continue 
c 

c define the elements 

c 

nmat=l 

ngen=0 

c 

do 30 j=l, jdmax-1 
do 31 k=l, kdmax-1 
do 32 i=l, idmax-1 

nelem=i+ (idmax-1) * (k-1) + ( idmax- 1) * (kdmax-1) * ( j-1) 

nl=i+idmax* (k-1) +idmax*kdmax* (j-1) 

n2=nl+l 

n3=n2+idmax 

n4-n3-l 

n5=nl+idmax*kdmax 
* n6-n5+l , 

n7==n6+idmax 
nS-nl-1 

c write (9,102) nelem, nmat , ngen, n5, n6, n7, n8, nl, n2, n3, n4 
c 102 format (1115) 

write (9, 102) nelem, nmat , n5, n 6, n7, n8, nl, n2, n3, n4 
102 format (i8, i5, 8i8) 

32 continue 
31 continue 
30 continue 
c 

nmat ^2 

c 

nelem=idmax-l+ (idmax-1) * (kdmax-1-1) 4- 
1 (idmax-1) * (kdmax-1) * ( jdmax-1 -1) 4-1 
nl=idmax4-idmax* (kdmax-1) 4-idmax*kdmax* ( jdmax-1) 4-1 
n2=nl+l 
n3=nl4-2 
n4=nl4-3 
n5=nl4-4 
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n6=nl+5 

n7=nl+6 

n8=nl+7 

c write ( 9, 10 4) nelem, nmat , ngen, nl, n2, n3, n4, n5, n6, n7 , n8 

c 104 format (lli5) 

write (9, 104) nelem, nmat, nl,n2, n3, n4, n5,n6, n7, n8 
104 format ( i8 , i5, 8i8 ) 
c 

c write (9, 110 ) 

c 110 format ( T 1 1 ll ? ,/, f 1 0 ',/,' 1 T ) 

c idmpl==idmax+l 

c idmp2 = idmax+2 

c write (9,111) idmp2 , idmpl, n4, n3, n2, nl 

c 111 f ormat ( 1 10 1 2 , ,2i5, 

cl /,' 1 0 1 , 4i5) 

c 

write (9, 110) 

110 format ( 1 1 1 ll T ,/, f 1 0\/, ? 

idmpl=idmax+l 

idmp2 = idmax+2 

write ( 9, 111) idmp2 , idmpl, n4, n3, n2, nl 

111 format ( T 1 1 2 f ,2i8, 

1 /,' l f ,4i8) 

c 

return 

end 

c 

subroutine rdavs 


c 

c**** subroutine rdavs 
c 


common /mesh/ il , i2,jl,j2,kl, k2, deltx, deity, deltz, 
1 xlow, ylow, zlow 


common/ cthdat/ velx (75, 75,75) , vely (75, 75, 75) , 

1 velz{75,75,75) , pres (75, 75, 75) , temp (75, 75, 75) , 

2 gama (75,75,75) , spen (75,75,75) , vmas (75,75,75), 

3 vphi (75,75,75) , cspd (75, 75,75) 


c 

common/cthxyz/xcth (100) , ycth (100) , zcth (100) 


1 T ) 


c 

c 


common /ve ling/ vfx (75,75,75) , vfy (75, 75, 75) , vfz (75, 75, 75) 
character*! chdum 


c 

open (2, file* 1 avsout 1 ) 
open (3, file* 1 avschk T ) 
c 

c calculate the nodal coordinates 

c 

do 20 i*il,i2 
vi=i 

xcth (i-il+1) =xlow+ (vi-1 . 0 ) *deltx-deltx 
20 continue 


do 21 j=jl, j2 
vj = j 

ycth ( j- j 1+1) =ylow+ (v j-1 . 0 ) *delty-delty 
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21 continue 
c 

do 22 k=kl,k2 
vk=k 

zcth (k-kl+1) =zlow+ (vk-1 . 0 ) *deltz-deltz 

22 continue 
c 

c read in the mesh data 

c 

do 400 idum=l , 27 
read ( 2 , 901) chdum 

901 format (al) 

400 continue 

c 

do 10 k=l,k2-kl+l 
do 11 j=l,j2-jl+l 
do 12 i—1, 12-il+l 

read (2, 902) velx(i, j,k) r vely (i, j,k) , velz (i, j, k) , 
1 pres (i, j, k) , temp(i, j,k) , gama (i, j, k) , spend, j,k) 
read (2, 902) vmas (i, j, k) , vphi (i, j, k) , 

1 cspd (i, j, k) 
c 

c unit conversion factors 

c 

velx ( i, j , k) =velx (i, j , k) /l . 0e+6 
vely (i, j , k) =vely (i, j, k) /l . 0e+6 
velz (i, j , k) =velz (i, j , k) /I . Oe+6 
pres (i, j , k) =pres (i, j , k) / I . (Te+12 
c 

902 format <7ell. 3) 

12 continue 

11 continue 
10 continue 
c 

do 30 k=l,k2-kl+l 
do 31 j=l , j2- j 1+1 
do 32 i=l,i2-il+l 

c write (3,903) velx (i, j, k) , vely (i, j, k) , velz (i, j , k) , 

c 1 pres (i, j, k) , temp (i, j, k) , gama (i, j,k) , spen (i, j, k) 
c. write (3,903) vmas (i, j , k) , vphi (i, j , k) , cspd (i, j , k) , 
c 1 xcth (i) , ycth ( j ) , zcth (k) 
c 903 format (7el 1.3) 

32 continue 
31 continue 
30 continue 
c 

c calculate the nodal velocities 
c 

do 40 k=2,k2-kl+l 
do 41 j~2, j2- jl+1 
do 42 i=2,i2-il+l 
wta=vmas (i, j, k-1) 
wtb=vmas (i, j-1, k-1) 
wtc=vmas (i, j-1, k) 
wtd=vmas (i, j, k) 
va=velx (i, j , k-1) 
vb=velx (i, j-1, k-1) 
vc=velx (i, j-1, k) 
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vd=velx (i f j, k) 

if ( wta+wtb+wtc+wtd. gt . 0 . 0 ) then 
vfx (i, j , k) = (wta*va+wtb*vb+wtc*vc+wtd*vd) / 
1 (wta+wtb+wtc+wtd) 

else 

vfx (i, j , k) =0 . 0 
endif 
c 

c write ( 10 , 910 ) wta, wtb, wtc, wtd, va, vb, vc, vd 

c 910 f ormat ( 8el0 . 2) 
c 

wt a=vmas (i-1, j , k-1) 
wtb=vmas ( i, j , k-1) 
wtc=vmas (i, j, k) 
wtd=vmas (i-1, j , k) 
va=vely (i-1, j, k-1) 
vb=vely (i, j, k-1) 
vc^vely (i, j, k) 
vd=vely (i-1, j, k) 

if (wta+wtb+wtc+wtd. gt . 0 . 0 ) then 
vf y ( i , j , k) = (wta*va+wtb*vb+wtc*vc+wtd*vd) / 
1 (wta+wtb+wtc+wtd) 

else 

vfy (i, j, k) -0 . 0 
endif 
c 

c write (10, 911) wta, wtb, wtc, wtd, va, vb, vc, vd 

c 911 format (8el0 .2) 
c 

wta=vmas (i, j , k) 
wtb=vmas ( i, j-1, k) 
wtc=vmas (i-1, j-1 , k) 
wtd=vmas (i-1, j , k) 
va=velz (i, j , k) 
vb=velz (i, j-1, k) 
vc=velz (i-1, j-1, k) 
vd^velz (i-1, j, k) 
if (wta+wtb+wtc+wtd. gt .0.0) then 
vf z (i, j , k) = (wta*va+wtb*vb+wtc*vc+wtd*vd) / 
1 (wta+wtb+wtc+wtd) 

else 

vf z (i, j, k) =0 . 0 
endif 
c 

c write ( 10, 912 ) wta, wtb, wtc, wtd, va, vb, vc, vd 

c 912 format (8el0 . 2) 
c 

42 continue 
41 continue 
40 continue 
c 

return 

end 

c 

subroutine dcelem(ii, j j, kk, 11) 
c 

c**** subroutine dcelem 
c 
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common/atype/itype, jtype 

common/ waldat/ywall, tcut , tmin, dencut , 

1 idmax, jdmax, kdmax, xdmin, xdmax, ydmin, ydmax, zdmin, zdmax 
common/props/den, denref , rmass, vmu, vlamb, eta, zeta, phi, r j (3, 3) 
common/ nda t a /xl (3) ,x2 (3) ,x3 (3) , x4 (3) , 

1 vl (3) , v2 (3) , v3 (3) , v4 (3) , xc<3) 

dimension h (3, 3) , f (3, 3) , e {3, 3) , ep(3,3),p(3) , 

1 y(33) , yt (33, 101) 
dimension param (50 ) , time (101) 


c 

c 

c 

c 

c 

c 


external fen 

call eldata (ii, j j , kk, 11) 

rden=den/ denref 
if (rden .gt . 1000 . 0) go to 99 
if (rden . It . 0 . 01) go to 99 
if (phi . gt . 0 . 999) go to 99 

denblk= (1 . -phi) *den/denref 
if (denblk . It . 0 . 01) go to 99 
if (denblk . It . dencut ) go to 99 


c 

c 


call nddata (ii, jj,kk,ll) 

call jcalc (den, xl, x2, x3, x4, xc, r j ) 
call jcalc (xl, x2, x3, x4, xc, rj) 


c 

c 

c 558 
c 


c 

c 557 
556 
- 555 

c 
c 
c 


c 

clOOO 

11 

10 

c 


c 

clOOl 


call initlz (h, f , ep, p) 
write ( 10 , 558) rmass 

format (/, 9 inertia tensor ( mass - \ell.3,' )',/) 

do 555 10=1,3 

do 556 jo=l,3 

write (10, 557) r j (io, jo) 

format (e20 . 8) 

continue 

continue 

set the state vector 

do 10 i=l,3 
do 11 j=l, 3 

y (3* (i-1) + j) =- h (i, j ) 

y<3*(i-l)+j+9> - f ( i , j ) 

y (3* (i-1) + j+18) - ep (i, j) 

write (10,1000) h(i,j),f(i,j),ep(i,j) 

format (3el5 . 4) 

continue 

continue 

do 12 i=l, 3 
y ( i+27 ) = p ( i) 
y (i+30) =xc (i) 
write (10,1001) p(i),xc(i) 
format (2el5 . 4) 


12 continue 
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c 

c set up and call the integration routine 

c 

ido=l 
neq=33 
tol=l . Oe-6 
c 

do 20 i=l, 50 
param { i) =0 . 0 
20 continue 
c 

param ( 1) =1 . 0e-6 
c param (3) =*1 . 0e-3 

param(4) =100000 
c param < 9) =1 . 0e-l 

c param(10)=3 

param < 1 2 ) —2 
c 

do 32 i=l / neq 
yt (i, 1) =y (i) 

32 continue 
c 

c iend=2 

iend=l 
riend=iend 
c 

tarr=l . 0e+33 
vceny-p (2) /rmass 
avceny=abs (vceny) 
if (a vceny .gt .0.0) 

1 tarr=abs ( (ywall-y (32) ) /vceny) 
c if (itype . eq. 2 . and. tarr . gt . tout ) go to 99 

c if (itype. eq. 2 .and.tarr.lt .tmin) go to 99 

c if (tarr .gt . tcut) go to 99 

if (tarr .ge . tcut) go to 99 
if (tarr . It . tmin) go to 99 
c 

c * if(ii.n$.3) go to 99 

c if(jj.ne.38) go to 99 

c if(kk.ne.44) go to 99 

c 

delt=tarr/riend 

c if (itype. eq. 2) delt=tarr/riend 

c if (itype .eq. 1) delt=0.2 

c 

t=0 . 0 

time (1) =0 . 0 
c 

do 30 i=l,iend 

c 

ri=i 

c t= (ri-1 . 0) *delt 

tend=t+delt 
c 

c call ivprk (ido, neq, fen, t , tend, tol, param, y) 

c call ivpbs (ido, neq, fen, t, tend, tol, param, y) 

c call ivpag (ido, neq, fen, f cn j , aaa, t , tend, tol, param, y) 
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C 

c call ifdcrk (ido, neq, t , tend, tol, param, y) 

call ivdcrk (ido, neq, t, tend, tol,param, y) 
c 

c write (10, 101) tarr, t , tend, y (32) 

c 101 format (5el5 . 3) 
c 

do 31 j=l,neq 
y t ( j , i+l) =y ( j ) 

31 continue 
c 

t=tend 

time (i+l) =tend 
c 

30 continue 
c 

c write cth input data 

c 

if (itype.eq. 1) 

1 call cthout (ii, j j , kk, 11, y) 
if (itype . eq. 2) 

1 call dynout (ii, j j , kk, 11, y, tarr) 
c 

ido=3 

c call ivprk (ido, neq, fen, t , tend, tol, param, y) 

c call ivpbs (ido, neq, fen, t, tend, tol, param, y) 

c call ivpag (ido, neq, fen, fen j , aaa, t, tend, tol, param, y) 

c 

c call prtout (lend, time, yt) 

c 

go to 99 
c 

write (7, 501) ii, j j, kk, 11, yt (31, 1) , yt (32, 1) , y t (33,1) 
write(8,501) ii,jj,kk,ll, 

1 yt (31, iend+1) , yt (32, iend+1) , yt (33, iend+1) 

501 format (4i6, 3el5 . 4) 
c 

99 return 
end 
c 

subroutine dynout (ip, jp, kp, Ip, y, tarr) 
c 

c **** subroutine dynout 
c 

common/props /den, denref , rmass, vmu, vlamb, eta, zeta, phi, r j (3,3) 
common/ndata/xl (3) , x2 (3) , x3 (3) , x4 (3) , 

1 vl (3) , v2 (3) , v3 (3) , v4 (3) ,xc (3) 

common/waldat/ywall, tcut , tmin, dencut, 

1 idmax, jdmax, kdmax, xdmin, xdmax, ydmin, ydmax, zdmin, zdmax 
c 

dimension y(33) ,f(3,3) , xtl (3) ,xt2(3) , xt3 (3) , xt4 (3) , xt5 (3) , 

1 ptc (3) , id (5) , jd (5) , kd (5) 
c 

cl= (denref /den) ** (-1 . /3 . ) 

c 

do 10 i=l, 3 
do 11 j=l,3 
f (i, j) =y (3* (i-1) +j + 9) 
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c write(4,101) i,j,f(i,j) 

c 101 format (15, 15, el5 . 6) 

11 continue 
10 continue 

c 

call detcal (f , detf ) 
denout=denref /detf 
c 

c momentum enhancement 

c 

ridmax=idmax 
r jdmax= jdmax 
rkdmax=kdmax 

tmass^denref * ( (xdmax-xdmin) / (ridmax-1 . 0) ) * 

1 ( (ydmax-ydmin) / (r jdmax-1 . 0) ) * < (zdmax-zdmin) / (rkdmax-1 . 0) ) 

enhmom= ( (tmass+rmass) / mass) **0.5 
c enhmom=l . 0 

c 

do 12 i=l f 3 

ptc (i) =enhmom*y (27+i) /5 . 0 

12 continue 
c 

do 20 i=l , 3 
xt 1 (i) =0 . 0 
xt2 (i) =0 . 0 
xt3 (i) =0 . 0 
xt 4 (i) =0 . 0 
xt5 (i) =0 . 0 

20 continue 
c 

do 23 i-1,3 
do 21 j-1,3 

xtl (i) -xtl (i) +cl*f (i, j ) * (xl ( j) -xc ( j) ) 
xt2 (i) ~xt2 (i) +cl*f (i, j) * (x2 (j) -xc ( j) ) 
xt3 (i) -xt3 (i) +cl*f (i, j) * (x3 < j) -xc ( j) ) 
xt4 (i) =xt 4 ( i) +cl*f ( i, j) * (x4 (j) -xc ( j) ) 

21 continue 
23 continue 

c 

do 22 i«l,3 
xtl(i)«xtl(i)+y(30+i) 
xt2 (i) -xt2 (i) +y (30+i) 
xt3 (i) =xt3 (i) +y (30+i) 
xt 4 (i) »xt4 (i) +y (30+i) 
xt5 (i) =y (30+i) 

22 continue 
c 

c calculate the element numbers 

c 

c nelmax- (idmax-1) * (jdmax-1) * (kdmax-1) 

c 

ridmax=idmax 

delx= (xdmax-xdmin) / (ridmax-1 . 0 ) 
c rjdmax= jdmax 

c dely= (ydmax-ydmin) / (r jdmax-1 . 0 ) 

rkdmax=kdmax 

delz- (zdmax-zdmin) / ( rkdmax-1. 0) 
c 
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id(l) 
kd(l) 
id (2) 
kd(2) 
id (3) 
kd<3) 
id ( 4) 
kd (4) 
id ( 5) 
kd (5) 


=int ( 
=int ( 
=int ( 
=int ( 
=int ( 
=int { 
=int ( 
=int ( 
=int ( 
=int ( 


(xtl(l) 
(xtl (3) 
(xt2 (1) 
(xt2 (3) 
(xt3 (1) 
(xt3 (3) 
(xt 4 ( 1 ) 
(xt4 (3) 
{xt 5 (1) 
(xt5 (3) 


■xdmin) 

-zdmin) 

■xdmin) 

■zdmin) 

■xdmin) 

■zdmin) 

•xdmin) 

•zdmin) 

•xdmin) 

■zdmin) 


/delx) +1 
/delz) +1 
/delx) -f-1 
/delz) +1 
/ delx) +1 
/delz) tl 
/delx) +1 
/delz) +1 
/delx) -HI 
/delz) -HI 


c 

c 

c 

c 


write(9,501) ip,jp,kp,lp 
501 format ( 4i5 ) 


do 303 ic=l,5 
ia=id (ic) 

if (ia .gt . idmax-1) go to 303 
if(ia.lt.l) go to 303 
ka=kd (ic) 

if (ka . gt . kdmax-1) go to 303 
if(ka.lt.l) go to 303 
nelem=ia-H (idmax-1) * (ka-1) 

write (9, 503) nelem,ptc (1) , ptc (2) , ptc (3) ,tarr 
write (61,503) nelem,ptc (1) , ptc (2) ,ptc (3) , tarr 
format (i5, 4el0 . 3) 
format (i8, 4el0 . 3) 
continue 

go to 301 


503 

503 

303 


c 

c 


write (9, 506) xtl (1) , xtl (2) f xtl (3) 

506 format (3e20 . 8) 

write (9, 507) xt2 (1) ,xt2 (2) ,xt2 (3) 

507 format (3e20 . 8 ) 

write (9, 50 8) xt3 (1) ,xt3 (2) ,xt3 (3) 

508 format (3e20 . 8) 

write (9, 598) xt4(l),xt4(2),xt4(3) 

598 format (3e20 . 8) 

write (9, 597) xt5 (1) , xt5 (2) , xt5 (3) 

597 format (3e20 . 8) 

301 dummy = 1 . 0 

return 

end 

subroutine cthout (ip, jp, kp, lp, y) 
r *** subroutine cthout 

common/props/den, denref, rmass, vmu, vlamb, eta, zeta,phi, r j (3,3) 
common/ ndata/xl (3) ,x2 (3) , x3 (3) , x4 (3) , 

1 vl (3) , v2 (3) , v3 (3) , v4 (3) , xc (3) 

dimension y (33) , f (3, 3) , xtl (3) , xt2 (3) ,xt3 (3) ,xt4 (3) , vtc (3) 

nummat=l 

numsub=49 



c 
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cl= (denref /den) ** (1 . /3 . ) 
c2= ( 1 . -phi) ** (1 . /3 . ) 
c3=c2/cl 
c 

do 10 i=l, 3 
do 11 j=l , 3 
f (i, j)»y(3*(i-l)+j + 9) 
c write ( 4 , 101) i,j,f(i,j) 

c 101 f ormat (15, i5, el5 . 6) 

11 continue 
10 continue 

c 

call detcal (f , detf ) 
denout=denref /detf 
c 

do 12 i=l, 3 

vtc(i)=1.0e+6*y(27+i) /rmass 

12 continue 
c 

do 20 i=l , 3 
xtl (i) =0 . 0 
xt2 (i) =0 . 0 
xt3 (i) =0 . 0 
xt4 (i) =0 . 0 

c write (4, 102) xl (i) , x2 (i) , x4 (i) , x4 (i) , xc (i) 

c 102 format (5el5 . 6) 

20 continue 
c 

do 23 i=l, 3 
do 21 j-1,3 

xtl (i) »xtl (i) +c3*f (i, j) * (xl ( j) -xc { j) ) 
xt2 (i) =xt2 (i) +c3*f (i, j) * (x2 (j) -xc ( j) ) 
xt3 (i) =xt3 (i) +c3*f (i r j) * (x3 ( j) -xc ( j) ) 
xt4 (i) =xt 4 (i) +c3*f (i f j) * (x4 ( j) -xc ( j) ) 

21 continue 
23 continue 

c 

c go to 301 
c 

do 22 i-1,3 
xtl (i) =xtl (i) +y (30 + i) 
xt2 (i) =xt2 (i) +y (30+i) 
xt3 (i) =xt3 (i) +y (30+i) 
xt4 (i) =xt4 (i) +y (30+i) 

22 continue 
c 

c 301 dummy= 1 . 0 
c 

write <4, 501) ip, jp, kp, lp 

501 format (5x, 'package ' , 13, ' , ' , 13, ' , ' , i3, ' , ' , i3) 
write (4, 502) nummat , numsub 

502 format (7x, 'material i3, /7x, 1 numsub 
write (4, 532) denout 

532 format (7x, 'density » r ,e20.8) 
write (4,503) vtc (1) 

503 f ormat (7x, ' xvel - ',e20.8) 
write (4,504) vtc (2) 


i3> 
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504 format (7x, 'yvel * T / e20.8) 
write(4,594) vtc(3) 

594 format (7x, ' zvel - *,e20.8) 
write ( 4 , 505) 

505 format (7x, ' insert pyramid 1 ) 
wr i t e ( 4 , 5 0 6 ) xt 1 ( 1 ) , xt 1 ( 2 ) , xt 1 ( 3 ) 

506 format ( 9x, 1 point = ',el5.6,', ',el5.6, 
write (4,507) xt2 ( 1) , xt2 (2 ) , xt2 ( 3) 

507 format ( 9x, 'point * ',615.6,*, T ,el5.6, 

write (4, 508) xt3 (1) , xt3 (2) , xt3 (3) 

508 format ( 9x, 'point = ',el5.6,', f ,el5.6, 

write (4, 598) xt 4 ( 1) , xt4 (2 ) , xt 4 (3) 

598 format (9x, 'vertex = ',el5.6,', ',el5.6, 

write (4,509) 

509 format (7x, 'endinsert ' ) 
write (4,591) 

591 format ( 5x, 'endpackage ' ) 
c 

return 
end 
c 

subroutine prtout (iend, time, yt) 
c 

c **** subroutine prtout 
c 

dimension yt(33,101) 
c 

dimension time (101) , out (18,101) ,£(3,3) , e (3, 3) , ep (3, 3) , 
1 c (3, 3) 
c 

c print output 

c 

c open ( 9, f ile- *pltdat/dc3d. out ' ) 

c 

open (11, f ile= 'pltdat/hll .pit 1 ) 
open (12, f ile= 'pltdat/h!2 .pit ' ) 
open (13, f ile= 'pltdat/h!3 .pit 1 ) 
open (14, f ile= 'pltdat/h21 .pit ' ) 

* open (15, f ile= *pltdat/h22 .pit 1 ) 
open (16, f ile~ 'pltdat/h23 .pit ' ) 
open (17, f ile= 'pltdat/h31 .pit ' ) 
open (18, f He- 'pltdat/h32 .pit ' ) 
open (19, f ile~ *pltdat/h33 .pit ' ) 
c 

open (20, f ile= 'pltdat/f 11 .pit * ) 
open (21, f ile= 'pltdat/f 12 .pit ' ) 
open (22, f ile= ' pltdat/f 13 .pit ' ) 
open (23, f ile= 'pltdat/f21 .pit ' ) 
open (24, f ile= ' pltdat/f 22 .pit ' ) 
open (25, f ile= 'pltdat/f 23 .pit ' ) 
open (26, f ile= 1 pltdat/f 3 1 .pit ' ) 
open (27, file= 'pltdat/f 32 .pit ' ) 
open (28, file= 'pltdat/f 33 .pit ' ) 
c 

open (2 9, file 35 'pltdat/epll .pit ' ) 
open (30, f ile= 'pltdat /epl2 .pit ' ) 
open (31, f ile= 'pltdat/epl3 .pit ' ) 
open (32, f ile= 'pltdat/ep21 .pit ' ) 


', ' , el5 . 6) 
', 1 , el5 . 6) 

T ,el5.6) 
', \el5.6) 
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open (33, file- 'pltdat/ep22 .pit 1 ) 
open (34, f ile= 'pltdat/ep23 .pit 1 ) 
open (35, f ile= ' pltdat /ep31 .pit T ) 
open (36, f ile= 'pltdat/ep32 .pit 1 ) 
open (37, f ile= 'pltdat/ep33 .pit 1 ) 
c 

open (38, f ile= 1 pltdat /pi .pit 1 ) 
open (3 9, f ile= 1 pltdat /p2 .pit ' ) 
open (40, f ile^ 1 pltdat/p3 .pit 1 ) 
c 

open (41, f ile= f pltdat /xcl .pit * ) 
open (42, f ile= T pltdat /xc2 .pit 1 ) 
open (43, f ile= r pltdat /xc3 .pit 1 ) 
c 

open (44, f ile= 'pltdat /detf .pit ' ) 
open (45, file= 'pltdat /ile .pit 1 ) 
open (4 6, f ile= 'pltdat/ j2e .pit ' ) 
open ( 47 , f ile= ' pltdat /ilep .pit 1 ) 
open (48, f ile= 'pltdat/ j2ep.plt * ) 
open (49, f ile= 1 pltdat /detc .pit 1 ) 
open (50, f ile= ' pltdat /out 7 .pit 1 ) 
open (51, f ile= 1 pltdat /out 8 .pit 1 ) 
open (52, f ile=* 1 pltdat /out 9 .pit 1 ) 
c 

c print out the states 
c 

c do 80 i=l,iend+l 
c write (9, 111) (yt ( j, i) , j*l, 33) 
c 111 format ( llel5 . 6) 
c 80 continue 
c 

c calculate output variables 
c 

do 60 i-1,18 
do 61 j=l,iend+l 
out (i, j) =0 . 0 
61 continue 
60 continue 
c 

do 610 i=l, iend+1 
c 

do 611 j=l, 3 
do 612 k=l f 3 

f(j,k) -yt<3*(j-l)+k+9,i) 
ep( j,k)-yt (3*( j-l)+k+18,i) 

612 continue 
611 continue 
c 

call detcal (f , detval) 
call ecalc(f,e) 
c 

do 614 j-1,3 
do 615 k-1,3 
if(j.eq.k) then 
del jk=l . 0 
else 

del jk=0 . 0 
endif 
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c ( j , k) =del jk+2 . 0*e < j, k) 

615 continue 
614 continue 
c 

call detcal (c, detc) 
c 

out <1, i) =detval 

out (2,i) = (l./3.) * (e(l, l)+e<2,2)+e(3,3> ) 
out (4,i) = (l./3.) * ( ep (1,1) +ep (2, 2) +ep (3, 3) ) 
out ( 6, i) =detc 
c 

610 continue 
c 

do 52 j=l, 9 

jl-j 

j2= j+9 
j3= j + 18 
kl-j+10 
k2- j+1 9 
k3= j+28 

do 51 i=l,iend+l 
write (kl, 112) time (i) , yt ( jl, i) 
write (k2 , 112) time (i) , yt ( j2, i) 
write (k3, 112) time (i) , yt ( j3, i) 

112 format (2e20 . 8) 

51 continue 

52 continue 
c 

do 53 j*l,3 
j 1= j+27 
j2= j+30 
kl= j+37 
k2= j+40 

do 54 i=l f iend+l 

write (kl, 102) time (i) , yt ( jl, i) 

write (k2, 102 ) time (i) , yt ( j2, i) 

102 f ormat (2e20 . 8 ) 

54 continue 

53 continue 
c_ 

do 55 j=l, 9 
j j=j+43 

do 56 i=l r iendtl 

write { j j, 113) time (i) , out ( j, i) 

113 format (2e20 . 8) 

56 continue 

55 continue 
c 

c do 50 i”l f iend+l 

c outl (i) ~yt (19, i) +yt (23, i) +yt (27, i) 
c out 2 (i) =abs (yt (19, i) -yt (23, i) ) 
c out3 (i) =abs (yt (23, i) -yt (27, i) ) 

c out 4 (i) =abs (yt (19, i) -yt (27 , i) ) 

c write ( 11, 101) time (i) , outl (i) , out2 (i) , out3 (i) , out 4 (i) 

c 101 format (5el5 . 6) 
c 50 continue 
c 


return 
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end 

c 

subroutine detcal (a, deta) 
c 

c**** subroutine detcal 
c 

dimension a (3, 3) 
c 

deta-a < 1, 1) * (a (2,2) *a (3, 3) -a (3,2) *a (2, 3) ) - 

1 a(l,2)*<a(2,l)*a(3,3)-a(2,3)*a<3,l) ) + 

2 a(l,3> *<a(2,l) *a<3,2) -a(3, 1) *a<2,2) ) 
c 

return 

end 

c 

subroutine eldata (ii, j j , kk, 11) 
c 

c**** subroutine eldata 
c 

common /props/ den , denref , mass, vmu, vlamb, eta, zeta, phi, r j (3, 3) 
c 

common /mesh/ il, i2, jl, j2, kl, k2 , del tx, deity, de It z, 

1 xlow, ylow, zlow 
c 

common/ cthdat/ ve lx (75, 75,75) , vely(75,75,75) , 

1 velz (75, 75, 75) , pres (75, 75, 75) , temp (75, 75, 75) , 

2 gama (75,75, 75), spen (7 5,75,75) , vmas (75,75,75) , 

3 vphi (75,75,75), cspd (75, 75,75) 
c 

c rmass=vmas (ii, j j , kk) 

c 

den=gama (ii, j j , kk) 
phi~vphi (ii, j j , kk) 
c 

c denref =2. 7 

c dmge=0 . 0 

c vmu * ( 1 . -dmge) *1 . 0 

c vlamb= ( 1 . -dmge) *1 . 0 

c eta=1.0 

c 

c write (10,100) velx (ii, j j , kk) , vely (ii, j j , kk) , velz (ii, j j , kk) , 

c 1 pres (ii, j j, kk) , temp (ii, j j , kk) , gama (ii, j j, kk) , spen (ii, j j , kk) , 
c 2 vmas (ii, j j , kk) , vphi (ii, j j , kk) , cspd(ii, j j, kk) 
c 100 format </7ell . 3, /3ell . 3, /) 
c 

return 

end 

c 

subroutine nddata (ii, j j , kk, 11) 
c 

c**** subroutine nddata 
c 

common/ndata/xl (3) ,x2(3) ,x3(3) ,x4(3), 

1 vl(3),v2(3),v3(3),v4(3),xc(3) 

c 

common/mesh/il, i2, jl, j2, kl, k2, deltx, deity, deltz, 

1 xlow, ylow, zlow 
c 
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common/ cthdat/ velx (75,75, 75) ,vely(75,75,75) , 

1 velz (75,75,75) ,pr es (75, 75, 75) , temp (75, 75, 75) , 

2 gama (75,75,75), spen (75,75,75), vmas (75,75,75) , 

3 vphi (75,75,75) , cspd (75, 75, 75) 
c 

common /cthxyz/xcth (100) , ycth (100) , zcth (100) 
c 

common /ve ling/ vf x (75, 75, 75 ) , vf y (75,75,75) ,vfz(75,75,75) 
c 

c set cell coordinates 

c 

xnl=xcth ( ii) 
xn2=xcth (ii+1) 
xn3=xcth (ii+1) 
xn4=xcth (ii) 
xn5-xcth (ii) 
xn6-xcth (ii+1) 
xn7=xcth (ii+1) 
xn8=xcth (ii) 
c 

ynl=ycth ( j j) 
yn2= s ycth ( j j) 
yn3=ycth ( j j + 1) 
yn4= s ycth ( j j+1) 
yn5=ycth ( j j ) 
yn6=ycth ( j j ) 
yn7-ycth ( j j+1) 
ynS^ycth ( j j+1) 
c 

znl~zcth (kk) 
zn2=zcth (kk) 
zn3=zcth (kk) 
zn4“zcth (kk) 
zn5=zcth (kk+1) 
znS^zcth (kk+1) 
zn7«zcth (kk+1) 
zn8=zcth (kk+1) 
c 

c set element coordinates 
- c 

go to (301,302,303,304,305) 11 
c 

301 xl (1) =xnl 
xl (2) =ynl 
xl (3) =znl 
x2 ( 1) ~xn2 
x2 (2)=yn2 
x2 (3) *zn2 
x3 (1) *xn4 
x3 (2) -yn4 
x3 (3) =zn4 
x4 ( 1) =xn5 
x4 (2) =yn5 
x4 (3) =zn5 
go to 306 

c 

302 xl(l)=xn3 
xl (2) =yn3 
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xl (3) =zn3 
x2 (1) =xn4 
x2 (2 ) =yn4 
x2 (3) = zn4 
x3 (1) =xn2 
x3 (2 ) =yn2 
x3 (3) =zn2 
x4 (1) =xn7 
x4 (2 ) =yn7 
x4 (3) =zn7 
go to 306 
c 

303 xl { 1 ) -xn4 
xl (2) =yn4 
xl (3) =zn4 
x2 ( 1) =xn2 
x2 <2)=yn2 
x2 (3) =zn2 
x3 (1) =xn7 
x3 (2) *=yn7 
x3 (3) =zn7 
x4 (1) =xn5 
x4 (2 ) =yn5 
x4 (3) =zn5 
go to 306 

c 

304 xl ( 1) ~xn8 
xl (2) =yn8 
xl (3) =zn8 
x2 ( 1) =xn5 
x2 (2 ) =yn5 
x2 (3) “zn5 
x3(l)«xn7 
x3 (2) =yn7 
x3 (3) =*zn7 
x4 ( 1) =xn4 
x4 (2) =yn4 
x4 (3) =zn4 

^ go to 306 
c 

305 xl (1) =xn6 
xl (2) =yn6 
xl (3) =zn6 
x2 (l)«xn7 
x2 (2 ) -yn7 
x2 (3) *zn7 
x3 ( 1) =xn5 
x3 (2) =yn5 
x3 (3) =zn5 
x4 (l)-xn2 
x4 (2 ) -yn2 
x4 (3) =zn2 
go to 306 

c 

306 dummy=l . 0 
c 

c set cell velocities 

c 
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vxnl=vf x ( ii, j j , kk) 
vxn2=vfx (ii+1, j j, kk) 
vxn3=vf x (ii+1, j j+1, kk) 
vxn4^vfx (ii, j j+1, kk) 
vxn5=vfx (ii, j j, kk+1) 
vxn6-vf x (ii+1, j j,kk+l) 
vxn7=vfx(ii, j j + 1, kk+1) 
vxn8=vf x (ii+1 , j j + 1 , kk+1) 
c 

vynl=vfy (ii, j j , kk) 
vyn2=vfy ( ii+1 , j j , kk) 
vyn3=vfy (ii+1, jj+l,kk) 
vyn4=vfy (ii, j j + 1 , kk) 
vyn5=vfy (ii, j j , kk+1) 
vyn6=vfy (ii+1, j j, kk+1) 
vyn7=vfy (ii, j j + 1, kk+1) 
vyn8=vf y ( ii+1, j j+1, kk+1) 

C 

vznl=vfz (ii, j j , kk) 
vzn2=vf z (ii + 1, j j , kk) 
vzn3-vfz (ii+1, j j+l,kk) 
vzn4=vf z ( ii, j j + 1 , kk) 
vzn5=vf z (ii, j j, kk+1) 
vzn6~vfz (ii+1, j j, kk+1) 
vzn7=vf z (ii, j j + 1, kk+1) 
vzn8=vfz (ii+1, j j+1, kk+1) 
c 

c set nodal velocities 

c 

go to (331,332,333,334,335) 11 

c 

331 vl ( 1 ) =vxnl 
vl (2 ) ^vynl 
vl (3) “vznl 
v2 (1) =vxn2 
v2 (2) =vyn2 
v2 (3) =vzn2 
v3 (1) =vxn4 
v3 (2) =vyn4 
v3 ( 3 ) =vzn4 
v4 ( 1 ) =vxn5 
v4 (2 ) =vyn5 
v4 (3) =vzn5 
go to 336 

c 

332 vl(l)=vxn3 
vl (2 ) =vyn3 
vl (3) =vzn3 
v2 ( 1) =vxn4 
v2 (2 ) =vyn4 
v2 (3) =vzn4 
v3 ( 1) =vxn2 
v3 (2) =vyn2 
v3 (3) =vzn2 
v4 (1) =vxn7 
v4 (2) =vyn7 
v4 (3) =vzn7 
go to 336 
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c 

333 vl ( 1) =vxn4 
vl (2) =vyn4 
vl (3) =vzn4 
v2 (1) =vxn2 
v2 (2) =vyn2 
v2 (3) =vzn2 
v3 (1) =vxn7 
v3 (2) =vyn7 
v3 (3) =vzn7 
v4 (1) =vxn5 
v4 (2) =vyn5 
v4 (3) =vzn5 
go to 336 

c 

334 vl ( 1 ) =vxn8 
vl (2 ) =vyn8 
vl (3) =vzn8 
v2 (1) =vxn5 
v2 (2 ) =vyn5 
v2 (3) =vzn5 
v3 ( 1) =vxn7 
v3 (2 ) =vyn7 
v3 ( 3) =vzn7 
v4 (1) =vxn4 
v4 (2 ) =vyn4 
v4 (3) =vzn4 
go to 336 

c 

335 vl (1) =vxn6 
vl ( 2 ) =*vyn6 
vl (3) =vzn6 
v2 ( 1 ) -vxn7 
v2 (2 ) =vyn7 
v2 ( 3 ) =vzn7 
v3 ( 1) *vxn5 
v3 (2) =vyn5 
v3 (3) =vzn5 
v4 ( 1) =vxn2 
v4 (2 ) =vyn2 
v4 ( 3 ) =vzn2 
go to 336 

c 

33 6 dummy =2 . 0 
c 

c set initial coordinates 

c 

go to 337 

xl (1) =0 . 0 
xl (2 ) =0 . 0 
xl (3) —0 . 0 
x2 (1)-1.0 
x2 (2 ) =0 . 0 
x2 (3) -0 . 0 
x3 ( 1) =0 . 5 
x3 (2) =sqrt (0 . 75) 
x3 (3) =0 . 0 


c 
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x4 (1) =0 . 5 

x4 (2) =sqrt (0 . 75) / 3 . 
x4 (3) =1 . 0 
c 

c set initial velocities 

c 

vl <1)=1.0 
vl (2 ) =0 . 0 
vl (3) =0 . 0 
v2 ( 1) =0 . 0 
v2 (2)~0.0 
v2 (3) =0 . 0 
v3 ( 1) =0 . 0 
v3 (2 ) =0 . 0 
v3 (3) =0 . 0 
v4 ( 1 ) =0 . 0 
v4 (2) =0 . 0 
v4 ( 3 ) =0 . 0 
c 

337 dummy =3 . 0 
c 

do 350 ip=l,3 

c write ( 10 , 3 40 ) xl (ip) , x2 (ip) , x3 (ip) , x4 (ip) 

c 340 f ormat (4el5 .*3) 

350 continue 
c write (10, 341) 
c 341 format { 1 1 ) 

do 352 ip=l,3 

c write (10 , 342) vl (ip) , v2 (ip) ,v3 (ip) , v4 (ip) 

c 342 format ( 4el5 . 3) 

352 continue 
c 

c 337 dummy =3 . 0 
c 

return 

end 

c 

subroutine initlz (h, f , ep, p) 
c 

c**** subrouitne initlz 
c 

common/props /den, denref , mass, vmu, vlamb, eta, zeta,phi, r j (3,3) 
common/ ndata/xl (3) , x2 (3) , x3 (3) , x4 (3) , 

1 vl (3) , v2 (3) , v3 (3) , v4 (3) ,xc(3) 

c 
c 

dimension h(3,3),f(3,3),ep(3,3),p<3) 
c 

c set initial conditions 
c 

f (1,1)=1.0 
f <l,2)-0.0 
f (1,3)=0.0 
f (2,1)=0.0 
f (2,2)=1.0 
f (2,3)=0.0 
f (3, 1)=0.0 
f (3,2)-0.0 
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C 


c 


c 

c 

c 

c 

c 

c 

c 

c 

c 103 
c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

-c 

c 

c 

c 

c 

c 

c 


f (3,3)=1.0 
detf=denref /den 

f (l,l)«detf**(l./3.) 
f (2,2)-detf**(l./3.) 
f <3,3>-detf**(l./3.) 

detf-f (1,1) *(f (2,2) *f (3,3)-f (3,2) *f (2,3) )- 

1 f (1,2) *(f <2, l)*f (3,3)-f (2,3) *f (3, 1) ) + 

2 f (1,3)* (f (2,1) *f (3,2) -f (3, l)*f (2,2)) 

den=denref /detf 

write (10, 103) den 
format (e20 . 8 ) 

ep (1, 1) =0 . 0 
ep { 1, 2) =0 . 0 
ep (1, 3) =0 . 0 
ep (2 , 1) =0 . 0 
ep <2 , 2 ) =0 . 0 
ep (2 , 3) -0 . 0 
ep (3, 1) =0 . 0 
ep (3, 2) =0 . 0 
ep (3, 3) =0 . 0 

h (1, 1) -0 . 0 
h (1, 2) =0 . 0 
h ( 1 , 3) =0 . 0 
h (2, 1) =0 . 0 
h (2 , 2 ) =0 . 0 
h (2, 3) *0 . 0 
h (3, 1) =0 . 0 
h (3, 2) *0 .0 
h (3, 3) =0 . 0 

p(l)“0.0 
p (2 ) =0 . 0 

p (3) -0 . 0 

call ivcalc(h,p) 

xc ( 1 ) =0 . 5 

xc (2) -sqrt (0 . 75) /3 . 
xc(3)=l. 0/4.0 

return 

end 


c 

subroutine f cn (neq, t , y, yprime) 
c 

c**** subroutine fen 
c 

common/props/den, denref , rmass, vmu, vlamb, eta, zeta, phi, r j (3,3) 
coiranon/ndata/xl (3) , x2 (3) , x3 (3) , x4 (3) , 

1 vl (3) ,v2 (3) ,v3 (3) ,v4 (3) , xc ( 3 ) 

c 


G-25 



dimension y (33) , yprime (33) 
c 

dimension h(3,3) r f (3 f 3),fi(3,3),e(3,3),ep(3 f 3) f p(3), 

1 rji(3,3),rm<3,3,3,3) ,rk(3,3) ,rkp(3,3) ,c<3,3) ,ci(3,3) , 

2 vg ( 3 , 3 ) , dd ( 3 , 3 ) , rn ( 3 , 3 , 3 , 3 ) , vkv (3,3) , xct ( 3 ) 
c 

dimension hd(3, 3) , fd<3, 3) ,ed<3, 3) ,epd<3, 3) ,pd(3) ,xcd(3) 
c 

c define variables 
c 

do 10 i-1, 3 
do 11 j-1, 3 
h(irj) - y(3*(i-l)+j) 
f (i, j) - y (3* ( i-1 ) + j + 9) 

©P(if j)“ y (3* (i-1) + j+18) 

11 continue 
10 continue 

c 

do 13 i-1, 3 
do 14 j-1, 3 
c (i, j) =0 . 0 
do 15 ia-1,3 

c (i, j) «c (i, j) +f (ia, i) *f (ia, j) 

15 continue 
14 continue 
13 continue 
c 

ninv-3 

ldc=3 

ldcinv=3 

c 

c call linrg (ninv, c, ldc, ci, ldcinv) 

if lag® 1 

call matv03 (if lag, ninv, c, ci) 
c 

do 12 i-1, 3 
P (i) * y (i+27) 
xct (i) = y (i+30) 

12 continue 
c 

call ecalc (f , e) 
c 

c write (10,104) 

c 104 format ( 1 subroutine fcn T ) 

c 

c call jcalc (den, xl, x2, x3, x4, r j ) 

c 

c write (10, 105) 

c 105 format ( 1 subroutine fen jcalc T ) 

c 

c do 201 ii-1,3 
c do 202 jj-1,3 
c write (10, 102) rj(ii,jj) 
c 102 format (e20 . 8) 
c 202 continue 
c 201 continue 
c 

call jinv (r j, r ji) 
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c 

c calculate transformer modulus 

c 

do 20 i=l,3 
do 21 j=l,3 
do 22 ia=l,3 
do 23 ib=l,3 
if(i.eq.ib) then 
delib=l . 0 
else 

delib=0 . 0 
endif 

if { j .eq. ib) then 
deljb-1.0 
else 

del jb^O . 0 
endif 

rm(i, j , ia, ib) =0 . 5*(delib*f (ia, j ) +f (ia, i) *del jb) 

23 continue 
22 continue 
21 continue 
20 continue 
c 

do 30 i=l,3 
do 31 j=l,3 
if(i.eq.j) then 
delij-1.0 
else 

deli j*0 . 0 
endif 
c 

c rk (i, j ) = (rmass/denref ) * (2 . *vmu* (e <i, j ) -ep (i, j) ) + 

c 1 deli j*vlamb* ( (e(l,l)+e(2,2)+e{3,3) ) - 
c 2 { ep { 1 , 1 ) +ep (2,2) +ep (3,3) ) ) ) 
c rkp (i, j) =-rk (i, j) 

c 

call detcal (f , rvl) 

rk ( i, j ) = (mass/ ( (l.-phi) *denref) ) *2.*vmu* ( (e(i, j) - 

1 < 1 . / 3 . ) *deli j* <e(l, 1) +e (2,2) +e (3, 3) ) )-ep(i, j) ) 

2 - (mass/ ( (1 . -phi) *denref ) ) * (vlamb+ (2 . /3 . ) *vmu) * (1 . /rvl-1 . ) 

3 *(l./rvl)*ci(i, j) 

rkp ( i , j ) — (mass/ ( (1 . -phi) *denref ) ) *2 . *vmu* ( (e (i, j) - 
1 (1./3.) *delij*(e(l,l)+e(2,2)+e(3,3> ) )-ep(i, j) ) 
c 

31 continue 
30 continue 
c 

do 50 i=l,3 
do 51 j-1,3 
fd(i, j>-0.0 
do 52 ia=l,3 

fd(i, j) =fd(i, j) +h (i, ia) *r ji (ia, j) 

52 continue 
51 continue 
50 continue 


c 


ninv=3 

ldf-3 
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ldf inv=3 
c 

c call linrg (ninv, f , Idf , f i, ldf inv) 

if lag=2 

call matv03 (if lag, ninv, f , f i) 
c 

do 300 i=l,3 
do 301 j=*l,3 
vg(i f j) =0 . 0 
do 302 ia=l,3 

vg ( if j ) =vg(i, j) +fd(i, ia) *fi (ia, j) 

302 continue 
301 continue 
300 continue 

c 

do 303 i-1,3 
do 304 j=l , 3 

dd(i, j) =0 . 5* (vg (i, j) +vg ( j , i) ) 

304 continue 

303 continue 
c 

c calculate transformer modulus 

c 

do 220 i-1,3 
do 221 j-1, 3 
do 222 ia»l f 3 
do 223 ib=l, 3 
if(i.eq.ib) then 
delib-1.0 
else 

delib-0 . 0 
endif 

if ( j .eq. ib) then 
del jb=l . 0 
else 

del jb=0 . 0 
endif 

rn (i, j f Lb, ia) =0 .5* (delib*fi (ia f j)+fi{ia, i) *deljb) 
223 continue 
^222 continue 
221 continue 
220 continue 
c 

c zeta=1.0 

call detcal (f , rvl) 
c 

do 231 i-1,3 
do 232 j-l f 3 
vkv (i r j ) =0 . 0 
do 233 ia-l f 3 
do 2 34 ib-l f 3 
if(ia.eq.ib) then 
delab=*l . 0 
else 

delab=0 . 0 
endif 

c vkv (i, j) =vkv (i, j) +zeta*rn (ia, ib, i, j) *dd(ia,ib) 

vkv (i, j ) =vkv {i, j ) +zeta* ( rmass*rvl/ ( (l.-phi) *denref ) ) 


i I 

u 
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1 

c 1 

c 2 

2 

3 

4 

5 

234 

233 

232 

231 

c 


43 

42 

41 

40 

c 


61 

60 

c 


70 

c 

c 

c 


81 

80 


82 


c 102 
202 
201 


*rn ( ia, ib, i, j) * 

delab* (1./3.) * (dd(l, 1) +dd<2, 2) +dd (3 f 3) ) 
delab* (1 . /3 . ) *aminl (0 . 0, dd(l, 1) +dd(2, 2) +dd(3 r 3) ) 
dd (ia, ib) 

* (0 .5* (dd(l, 1) **2+dd(l f 2) **2+dd(l,3) **2 + 
dd(2 r 1) **2+dd(2 f 2) **2+dd(2/ 3) **2 + 
dd ( 3 , 1 ) **2+dd<3/2) **2+dd(3 f 3) **2) ) **2 

continue 

continue 

continue 

continue 

do 40 1=1,3 
do 41 j=l , 3 
hd (i, j) =0 . 0 
do 42 ia=l,3 
do 43 ib=l,3 

hd (i, j) =hd(i, j) -rm(ia, ib f i, j ) *rk (ia f ib) 

continue 

continue 

hd(i, j ) =hd (i, j) -vkv{i, j) 

continue 

continue 

do 60 1=1/3 
do 61 j=l,3 

epd(i r j) =- (1. /eta) *rkp (i f j) 

continue 

continue 

do 70 1=1/3 
pd ( i) =0 . 0 

xcd (1) = (1 . /rmass) *p (i) 
continue 

set the vector yprime 

do 80 i=l/3 
do 81 j=l,3 

yprime (3* (i-1) +j) - hd(i/j) 

yprime (3* (1-1) +j + 9) = fd(i/j) 

yprime (3* (i-1) +j+18) = epd(i/j) 
continue 
continue 

do 82 1=1/3 
yprime (i +27) = pd(i) 
yprime (i+30) =xcd(i) 
continue 

do 201 ii=l f 3 
do 202 j j=l, 3 

write (10/ 102) f (11/ j j) f fd(ii/ j j) f h (11/ j j) , hd(ii f j j) 

format (5el5 . 3) 

continue 

continue 

return 
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end 

c 

subroutine ecalc(f,e) 
c 

c**** subroutine ecalc 
c 

dimens ion f(3,3),e(3,3),c(3,3) 
c 

c calculate e 

c 

do 10 i=l , 3 
do 11 j-1,3 
c 

c e (i, j) =0 . 0 

c (i, j ) =0 . 0 
if (i .eq. j) then 
deli j 38 ! . 0 
else 

deli j=0 . 0 
endif 
c 

do 12 ia=l, 3 

c (i f j) =c {i, j ) +f (ia, i) *f (ia, j) 

12 continue 
c 

e (i, j) =0 .50* (c (i, j) -deli j) 
c 

11 continue 
10 continue 
c 

return 

end 

c 

subroutine jinv(z,zi) 
c 

c **** subroutine jinv 
c 

c inverts the symmetric 3x3 matrix z 

c 

dimension z (3, 3) , zi (3, 3) 
c 

zl=z (1, 1) 
z2*z (2 , 2 ) 
z3-z (3, 3) 
z4=z (1, 2) 
z5~z (2, 3) 
z6=z (1,3) 
c 

det=zl* (z2*z3-z5*z5) -z4* (z4*z3-z5*z6)+z6* (z4*z5-z2*z6) 
c 

zi(l, 1) =+ (1 . /det) * (z2*z3-z5*z5) 
zi (2, 1) =- (1 . /det) * (z4*z3-z5*z6) 
zi<3, 1) =+ (1 . /det) * ( z4*z5-z2 *z 6) 
zi (1, 2) — (1 . /det) * (z4*z3-z5*z6) 
zi ( 2 , 2 )=+<!. /det) *(zl*z3-z6*z6) 
zi{3,2) — (1 . /det) * (zl*z5-z4*z6) 
zi <1, 3) =+ (1 . /det) * (z4*z5-z2*z6) 
zi (2, 3) (1 . /det) * < zl*z5-z4*z6) 


Ill ill Mini 
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zi (3, 3) =+ (1 . /det) *(zl*z2-z4*z4) 
c 

return 

end 

c 

subroutine jref(g,rjp) 
c 

c**** subroutine jref 
c 

c ommon /props/ den , den ref, rmass, vmu, vlamb, eta, zeta, phi, r j ( 3 , 3) 
common/rdata/xxl (3) , xx2 (3) , xx3 (3) , xx4 (3) , xxc (3) , vref 
c 

dimension g (3, 3) , r jp (3, 3) 
c 

detg=g (1, l)*(g(2,2)*g(3,3)-g(3,2)*g(2,3))- 

1 g<l,2) *(g(2,l) *g (3, 3) -g (2 , 3) *g (3, 1) ) + 

2 g(l,3) Mg (2,1) *g (3, 2) -g (3, 1) *g (2, 2) ) 
c 

c rmass= (1 . 0-phi) *den*vref *detg 

rmass=abs { (1 . 0-phi) *den*vref *detg) 
c 

c calculate the reference inertia tensor 

c 

hi = sqrt ( (xxl (1) -xxc (1) ) **2+ (xxl (2) - 
1 xxc (2) ) **2+ (xxl (3) -xxc (3) ) **2) 
h2 = sqrt ( <xx2 ( 1) -xxc (1) ) **2+ (xx2 (2) - 
1 xxc (2) ) **2+ (xx2 (3) -xxc (3) ) **2 ) 
h3 ** sqrt ( (xx3 ( 1) -xxc ( 1) ) **2+ (xx3 (2) - 
1 xxc (2) ) **2+ (xx3 (3) -xxc (3) ) **2) 
h4 = sqrt ( (xx4 (1) -xxc (1) ) **2+ (xx4 (2) - 
1 xxc (2) ) **2+ (xx4 (3) -xxc (3) ) **2) 
c 

epll - (xxl (1) -xxc(l) ) /hi 
epl2 = (xxl (2 ) -xxc (2 )) /hi 
epl3 = (xxl (3) -xxc (3) ) /hi 
c 

ep21 = (xx2 (1) -xxc(l) ) /h2 
ep22 = (xx2 (2 ) -xxc (2) ) /h2 
- ep23 = (xx2 (3) -xxc (3) ) /h2 
c 

ep31 = (xx3 (1) -xxc ( 1) ) /h3 
ep32 - (xx3 (2 ) -xxc (2) ) /h3 
ep33 = (xx3 (3) -xxc (3) ) /h3 
c 

ep41 * (xx4 (1) -xxc (1) ) /h4 
ep42 * (xx4 (2 ) -xxc (2 ) ) /h4 
ep43 = (xx4 (3) -xxc (3) ) /h4 
c 

r jpll= ( rmass /20 . ) * ( (hi* *2) *epll*epll+ (h2**2) *ep21*ep21 
1 +<h3**2) *ep31*ep31+ (h4**2) *ep41*ep41) 
r jpl2= ( rmass /20 . ) * ( (hl**2) *epll*epl2+ (h2**2 ) *ep21*ep22 
1 + (h3**2 ) *ep31*ep32+ (h4**2) *ep41*ep42) 
r jp!3= (rmass/20 . ) * ( (hi* *2) *epll*epl3+ <h2**2) *ep21*ep23 
1 + (h3**2 ) *ep31*ep33+ (h4**2 ) *ep41*ep43) 


r jp2 1= ( rmass/20 . ) * ( (hl**2) *epl2*epll+ (h2**2) *ep22*ep21 
1 + (h3**2) *ep32*ep31+ (h4**2) *ep42*ep41) 
r jp22= (rmass/20 . ) * ( (hl**2) *ep!2*ep!2+ (h2**2) *ep22*ep22 



1 + (h3**2 ) *ep32*ep32+ (h4**2 ) *ep42*ep42) 
r jp23= ( rmass/20 . ) * ( (hi* *2) *epl2*ep!3+ (h2**2) *ep22*ep23 
1 +(h3**2) *ep32*ep33+(h4**2) *ep42*ep43) 

r jp31= (rmass/20 . ) * ( (hi* *2) *epl3*epll+ (h2**2) *ep23*ep21 
1 + (h3**2 ) *ep33*ep31+ (h4**2) *ep43*ep41) 
r jp32= (rmass/20 . ) * ( (hi* *2) *ep!3*epl2+ (h2**2) *ep23*ep22 
1 + (h3**2) *ep33*ep32+ (h4**2) *ep43*ep42) 
rjp33=( rmass/2 0.) * ( (hi* *2) *epl3 *epl3-f- (h2**2 ) *ep2 3*ep2 3 
1 + (h3**2) *ep33*ep33+ (h4**2 ) *ep43*ep43) 

r jp (1, 1) =r jpll 
r jp (1, 2) =r jpl2 
r jp ( 1 , 3) -r jpl3 
r jp (2, 1) ~r jp21 
r jp (2,2) =r jp22 
rjp(2,3)=rjp23 
r jp (3, 1) =r jp31 
r jp (3,2) =r jp32 
r jp (3, 3) =r jp33 

return 

end 

subroutine jcalc (xl, x2, x3, x4 f xc, r j) 
subroutine jcalc (den, xl, x2, x3, x4, xc, r j) 

subroutine jcalc 

common/ rdat a /xxl (3) , xx2 (3) , xx3 (3) , xx4 (3) , xxc (3) , vref 

dimension xl (3) , x2 (3) , x3 (3) , x4 (3) , r j (3, 3) , 

1 r jp (3,3) , a (12, 12) ,b (12) , s (12) , ainv (12, 12) , 

2 g (3, 3) , ginv (3,3) , xc (3) , ami (12,12), am2 (12,12) 

set the reference tetrahedron 

xxl (1)=0.0 

xxl (2) -0 . 0 

xxl (3) -0 . 0 

xx2 (l)-l.O 

xx2 (2) -0.0 

xx2 (3) =0.0 

xx3 ( 1) =0 . 5 

xx3 ( 2 ) =sqrt (0.75) 

xx3 (3) “0 . 0 

xx4 (1) =0 . 5 

xx4 (2) =sqrt (0 .75) /3 . 

xx4 (3) =1 . 0 

xxc ( 1) =0 . 5 

xxc (2 ) =sqrt (0.75) /3 . 

xxc (3) =1 . 0/4 . 0 

vref=(l./3.) * (1./2.) *sqrt (0.75) 
set the right hand side vector 
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b(i) -xl(l) 
b(2) =xl (2) 
b (3) =xl (3) 
b (4) =x2 (1) 
b (5) =x2 (2) 
b(6) =x2 (3) 
b (7) =x3 ( 1) 
b (8 ) =x3 (2 ) 
b ( 9) =x3 (3) 
b ( 10 ) =x4 (1) 
b ( 11) =x4 (2) 
b ( 12 ) =x4 (3) 
c 

c set the coefficient matrix 

c 

a (1, 1) =xxl (1) -xxc (1) 
a (1,2) =xxl (2) -xxc (2 ) 
a (1,3) =*xxl (3) -xxc (3) 
a (1, 4) =0.0 
a ( 1, 5) =0.0 
a(l,6) =0.0 
a(l,7) =0.0 
a (1, 8 ) =0.0 
a (l f 9) =0.0 
a (1, 10) =1 . 0 
a (1, 11) =0 . 0 
a ( 1, 12 ) =0 . 0 
c 

a (2 , 1 ) =0.0 
a (2,2) =0.0 
a (2 , 3) =0.0 
a (2, 4) =xxl (1) -xxc (1) 
a < 2 , 5 ) =xxl ( 2 ) -xxc ( 2 ) 
a (2,6) =xxl (3) -xxc (3) 
a (2, 7) =0.0 
a (2, 8 ) =0.0 
a (2, 9) =0.0 
a (2, 10) =0 . 0 
a (2, 11) =1 . 0 
a (2 , 12 ) =0 . 0 
c 

a (3, 1) =0.0 
a (3, 2 ) =0.0 
a (3,3) =0.0 
a (3, 4) =0.0 
a (3, 5) =0.0 
a (3, 6) =0.0 
a ( 3 , 7 ) =xxl { 1 ) -xxc ( 1 ) 
a (3, 8) =xxl (2) -xxc (2 ) 
a ( 3 , 9 ) =xxl ( 3 ) -xxc ( 3 ) 
a (3, 10) =0 . 0 
a (3, 11) =0 . 0 
a (3, 12) =1 . 0 

a ( 4 , 1 ) =xx2 ( 1 ) -xxc ( 1 ) 
a ( 4 , 2 ) =xx2 ( 2 ) -xxc ( 2 ) 
a ( 4 , 3 ) =xx2 ( 3 ) -xxc ( 3 ) 
a (4, 4) =0.0 


< 2-2 
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} I i 


a ( 4, 5) =0.0 
a (4, 6) =0.0 
a ( 4 f 7) =0.0 
a (4 f 8) =0.0 
a (4 , 9) =0.0 
a (4, 10) =1 . 0 
a (4, 11) =0 . 0 
a (4, 12) =0 . 0 

a (5, 1) =0.0 
a (5, 2 ) =0.0 
a (5, 3) =0.0 
a ( 5 , 4 ) =xx2 ( 1 ) -xxc ( 1 ) 
a ( 5 , 5 ) =xx2 ( 2 ) -xxc ( 2 ) 
a ( 5 , 6 ) =xx2 ( 3 ) -xxc ( 3 ) 
a (5, 7 ) =0.0 
a<5,8) =0.0 
a (5, 9) =0.0 
a (5, 10) =0 . 0 
a (5, 11) =1 . 0 
a (5, 12) =0 . 0 

a ( 6, 1) =0.0 
a ( 6, 2 ) =0.0 - 
a(6,3) =0.0 
a ( 6, 4) =0.0 
a ( 6, 5) =0.0 
a ( 6, 6) =0.0 
a (6, 7) =xx2 (1) -xxc ( 1) 
a (6, 8) =xx2 (2 ) -xxc {2 ) 
a ( 6 , 9 ) =xx2 ( 3 ) -xxc ( 3 ) 
a ( 6, 10) =0 . 0 
a ( 6, 11) =0 . 0 
a (6, 12) =1 . 0 

a ( 7 f 1 ) =xx3 ( 1 ) -xxc ( 1 ) 
a(7,2) =xx3 (2) -xxc (2 ) 
a ( 7 , 3 ) =xx3 ( 3 ) -xxc ( 3 ) 
a (7, 4) =0.0 
a (7, 5) =0.0 
a (7, 6) =0.0 
a (7 , 7) =0.0 
a (7, 8 ) =0.0 
a (7, 9) =0.0 
a {7, 10) =1 . 0 
a (7, 11) =0 . 0 
a(7, 12) =0 . 0 

3(8,1) =0.0 
a ( 8, 2 ) =0.0 
a (8 , 3) =0.0 
a (8, 4) =xx3 (1) -xxc ( 1) 
a (8, 5) =xx3 (2 ) -xxc (2 ) 
a ( 8 , 6 ) =xx3 ( 3 ) -xxc ( 3 ) 
a(8,7) =0.0 
a(8,8) =0.0 
a <8 , 9) =0.0 
a (8 , 10) =0 . 0 


i ; 

u 


i til 


- t 



L 



tl 

S I 

1_.J 


a (8, 11)=1.0 


c 

a (8, 12) =0 . 0 

- - 


a (9, 1) =0.0 

— * 


a (9, 2 ) =0.0 

_ , 


a ( 9, 3) =0.0 
a (9, 4) =0.0 

= 


a ( 9, 5) =0.0 



a ( 9, 6) =0.0 

5. 

i 1 


a (9, 7) =xx3 (1) -xxc (1) 
a (9, 8) =xx3 (2) -xxc (2) 



a (9, 9) =xx3 (3) -xxc (3) 


c 

a (9, 10)=0 .0 
a ( 9, 11) =0 . 0 
a (9, 12) =1 . 0 



a (10,1) =xx4 (1) -xxc (1) 
a (10,2) =xx4 (2) -xxc (2) 

tj 


a(10,3) =xx4 (3) -xxc (3) 



a (10, 4) =0.0 

„ . 


a ( 10 , 5) =0.0 
a (10, 6) =0.0 

- - 


a (10, 7) =0.0 

w 


a (10, 8) =0.0 

i :f 


a (10, 9) =0.0 
a (10, 10)=1.0 

l;| 


a (10, 11) =0 . 0 


c 

a (10, 12) =0 . 0 

1 ; 


a (11, 1) =0.0 



a (11, 2) =0.0 

1 i 


a (11, 3) =0.0 
a (11, 4) =xx4 (1) -xxc (1) 



a (11, 5) =xx4 (2) -xxc (2) 

— 


a (11, 6) — xx4 (3) -xxc (3) 

- 


a (11, 7) =0.0 
a (11, 8) =0.0 



a (11, 9) =0.0 

y 


a (11, 10)=0.0 


- 

a (11, 11) =1 . 0 

i ( 


a (11, 12) =0 . 0 

— 

c 




a (12, 1) =0.0 



a (12, 2) =0.0 
a (12, 3) =0.0 



a (12, 4) =0.0 

ti 


a (12, 5) =0.0 
a (12, 6) =0.0 


a (12, 7) =xx4 (1) -xxc (1) 

y 


a(12,8) =xx4 (2) -xxc (2) 



a (12, 9) =xx4 (3) -xxc (3) 
a (12, 10) =0 . 0 
a (12, 11) =0 . 0 



a (12, 12) =1 . 0 
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call inverse solver 


c 

c 

c 


ninv=12 
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lda=12 

ldainv=12 

c 

c write (10,106) 

c 106 format ( ? sub jcalc before solver 1 ) 
c 

c do 108 iii=l, 12 
c do 109 j j j = l, 12 
c write ( 10 , 110 ) a(iii,jjj) 

c 110 format (e20 . 8) 
c 109 continue 
c 108 continue 
c 

c call linrg (ninv, a, Ida, ainv, ldainv) 

if lag=3 

call matvl2 (if lag, ninv, a, ainv) 
c call amodl(a,aml) 

c call mat inv (ninv, ami, am2) 

c call amod2 (am2 , ainv) 

c 

c write (10, 107) 

c 107 format ('sub jcalc after solver') 
c 

c calculate the solution vector 

c 

do 10 i-1,12 
s (i) =0 . 0 
do 11 j=l, 12 
s (i) =s (i) +ainv (i, j) *b ( j) 

11 continue 

c write(4, 876) i,b(i),s(i) 
c 876 format (lx, i5, 2el5 . 3) 

10 continue 
c 

c set the tensor g 

c 

g (1, 1) =s (1) 
g (1, 2) =s (2) 

* g (1, 3) =s (3) 
g (2, 1) =s (4) 
g (2 , 2) =s (5) 
g<2,3)-s(6) 
g (3, 1) =s (7) 
g (3, 2) =s (8) 
g (3, 3) =s ( 9) 
c 

xc (l)=s (10) 
xc (2)=s (11) 
xc (3) =s ( 12 ) 
c 

c invert the tensor g 

c 

c call inverse solver 

c 

ninv=3 

ldg=3 

ldginv^O 

c 
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c call linrg (ninv, g, ldg, ginv, ldginv) 

if lag=4 

call matv03 ( if lag, ninv, g, ginv) 
c 

call jref (g, r jp) 
c 

c calculate the inertia tensor 

c 

do 20 i=l,3 
do 21 j-1,3 
r j (i, j) =0 . 0 
do 22 ia=l,3 
do 23 ib=l, 3 

r j (i, j) =r j (i, j ) +ginv (ia, i) *rjp(ia,ib) *ginv(ib, j) 

23 continue 
22 continue 
21 continue 
20 continue 
c 

return 

end 

c 

c 

subroutine ivcalc(h,p) 
c 

c**** subroutine ivcalc 
c 

common /props /den, denref , rmass, vmu, vlarab, eta, zeta, phi, r j (3,3) 
common/ rdata/xxl (3) , xx2 (3) , xx3 (3) , xx4 (3) , xxc (3) , vref 
common/ ndata/xl (3) , x2 (3) , x3 (3) , x4 (3) , 

1 vl (3) , v2 (3) ,v3 (3) ,v4 (3) ,xc(3) 

c 

dimension a (12, 12) ,b (12) , s (12) , ainv (12, 12) , 

1 g(3,3) ,h(3,3) ,p(3) , ami (12,12) , am2 (12,12) 
c 

c set the right hand side vector 

c 

b (1) = v 1(1) 
b (2) =vl (2 ) 
b (3) =vl (3 ) 
b (4) =v2 ( 1 ) 
b (5) =v2 (2) 
b (6) =v2 (3) 
b (7) =v3 (1) 
b (8) =v3 (2) 
b ( 9) =v3 (3) 
b (10) =v4 (1) 
b(ll)=v4 (2) 
b ( 12 ) =v4 (3) 
c 

c set the coefficient matrix 

c 

a ( 1 , 1 ) =xl(l)-xc(l) 
a (1, 2) =xl (2) -xc (2) 
a ( 1, 3) =xl (3) -xc (3) 
a (1, 4) =0.0 
a ( 1, 5) =0.0 
a (1, 6) =0.0 
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a (1, 7) =0.0 
a ( 1 , 8 ) =0.0 
a ( 1 , 9) =0.0 
a (1, 10) =1 . 0 
a (l f 11) =0 . 0 
a ( 1 , 12 ) =0 . 0 

a (2, 1) =0.0 
a (2,2) =0.0 
a (2 , 3 ) =0.0 
a (2 , 4 ) =xl (1) ~xc (1) 
a (2, 5) =xl (2) -xc (2) 
a (2, 6) =xl (3) -xc (3) 
a (2,7) =0.0 
a < 2 , 8 ) =0.0 
a ( 2 , 9 ) =0.0 
a (2 , 10) =0 . 0 
a (2, 11) =1 . 0 
a (2, 12 ) =0 . 0 

a (3, 1) =0.0 
a (3, 2) =0.0 
a (3, 3) =0.0 
a (3, 4) =0.0 
a (3, 5) =0.0 
a (3, 6) =0.0 
a (3, 7 ) =xl ( 1 ) — xc ( 1 ) 
a(3,8) =xl (2) -xc (2) 
a (3, 9) =xl (3) -xc (3) 
a (3, 10) =0 . 0 
a (3, 11) =0 . 0 
a (3, 12) =1 . 0 

a (4, 1) =x2 (1) -xc (1) 
a {4, 2) =x2 (2) -xc (2) 
a (4, 3) =x2 (3) -xc (3) 
a (4, 4) =0.0 
a (4 f 5) =0.0 
a (4, 6) =0.0 
a (4,7) =0.0 
a (4,8) =0.0 
a ( 4 , 9) =0.0 
a (4, 10) =1 . 0 
a (4 , 11) =0 . 0 
a(4,12)-0.0 

a {5, 1) =0.0 
a(5,2) =0.0 
a (5, 3) =0.0 
a (5 r 4) =x2 (1) -xc (1) 
a (5, 5) =x2 (2) -xc (2) 
a ( 5 , 6 ) =x2 ( 3 ) -xc ( 3 ) 
a (5, 7) =0.0 
a (5, 8) =0.0 
a (5, 9) =0.0 
a(5, 10) =0 . 0 
a(5, 11)=1.0 
a (5, 12 ) =0 . 0 
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C 




a(6,l) -0.0 
a (6,2) =0.0 



a ( 6, 3) =0.0 



a ( 6, 4) =0.0 
a ( 6, 5 ) =0.0 
a ( 6, 6) =0.0 

— 


a ( 6, 7) =x2 (1) -xc (1) 
a ( 6, 8) =x2 (2) -xc (2) 
a ( 6, 9) =x2 (3) -xc (3) 

u 

C 

a (6, 10) =0 . 0 
a (6, 11)=0.0 
a ( 6, 12 ) =1 . 0 

S' i 


a (7, 1) =x3 (1) -xc (1) 



a (7, 2) =x3 (2) -xc (2) 
a(7,3) =x3(3)-xc(3) 
a (7, 4) =0.0 

- ■ 


a (7, 5) =0.0 



a (7, 6) =0.0 
a (7, 7) =0.0 

u 


a (7 , 8 ) =0.0 
a (7, 9) =0.0 
a (7, 10) =1 . 0 • 

- 


a (7, 11) =0 . 0 
a (7, 12) =0 . 0 

— 

C 

a (8, 1) =0.0 



a (8, 2) =0.0 

w ~ 1 


a (8, 3) =0.0 

1 1 


a (8, 4) =x3 (1) -xc (1) 
a (8, 5) =x3 (2) -xc (2) 



a ( 8 , 6 ) =x3 ( 3 ) -xc ( 3 ) 



a (8, 7) =0.0 



a (8, 8) =0.0 
a (8, 9) =0.0 

i ' : i 


a (8, 10) =0 . 0 



a (8, 11) =1 . 0 



a (8, 12) =0 . 0 


c 

a (9, 1) =0.0 



a (9, 2) =0.0 
a (9, 3) =0.0 

y 


a (9, 4) =0.0 
a (9, 5) =0.0 
a < 9, 6) =0.0 
a (9, 7) =x3 (1) -xc (1) 

i i 


a (9, 8) =x3 (2 ) -xc (2 ) 

*=3 


a (9, 9) =x3 (3) -xc (3) 



a (9, 10) =0 . 0 
a { 9, 11) =0 . 0 

' - 


a (9, 12 ) =1 . 0 

| J 

c 




a (10, 1) =x4 (1) -xc (1) 
a (10,2) =x4 (2) -xc (2) 

IrS 

u 


a (10, 3) =x4(3)-xc(3) 
a (10, 4) =0.0 
a (10, 5) =0.0 
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r ? 




a (10 f 6) =0.0 
a (10,7) -0.0 
a ( 10 , 8) =0.0 
a (10, 9) =0.0 
a (10, 10) *1 . 0 
a (10, 11) =0 . 0 
a (10, 12) =0 . 0 
c 

a (11, 1) =0.0 
a ( 1 1 , 2 ) =0.0 

a (11, 3) =0.0 
a (11, 4) =x4 (1) -xc (1) 

a ( 11 , 5) =x4 (2) -xc(2) 
a (11, 6) =x4(3) -xc (3) 
a ( 11 , 7 ) =0.0 

a ( 11, 8 ) =0.0 

a ( 1 1 , 9 ) =0.0 

a ( 11 , 10 ) =0 . 0 
a (11, 11) =1 . 0 
a (11, 12) =0 . 0 
c 

a ( 12 , 1) =0.0 
a ( 12, 2 ) =0.0 
a ( 12 , 3) =0.0 
a ( 12 , 4) =0.0 
a ( 12 , 5) =0.0 
a ( 12 , 6) =0.0 
a (12, 7) =x4 (1) -xc ( 1 ) 
a ( 12 , 8 ) =x4 (2) -xc (2 ) 
a (12, 9) =x4 ( 3 ) -xc ( 3 ) 
a (12, 10)=0 .0 
a (12, 11) =0 . 0 
a (12 , 12) -1 . 0 
c 

c call inverse solver 

c 

ninv=12 
lda=12 
ldainv=12 
- c 

c write ( 10 , 106) 

c 106 format ('sub jcalc before solver') 
c 

c do 108 iii=l, 12 

c do 109 jjj=l,12 
c write (10, 110) a(iii,jjj) 

c 110 format (e20 . 8) 
c 109 continue 
c 108 continue 
c 

c call linrg (ninv, a, Ida, ainv, Idainv) 

if lag=5 

call matvl2 (if lag, ninv, a, ainv) 
c call amodl(a,aml) 

c call mat inv (ninv, ami, am2 ) 

c call amod2 (am2, ainv) 

c 
c 


write (10,107) 
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c 107 format ('sub jcalc after solver 1 ) 
c 

c calculate the solution vector 

c 

do 10 i-1 , 12 
s (i) =0 . 0 
do 11 j-1,12 
s (i) =s (i) +ainv (i, j) *b( j) 

11 continue 
10 continue 
c 

c set the tensor g 

c 

const 35 (denref /den) ** (1 . /3 . ) 
c 

g (1, 1) =s (1) *const 
g ( 1, 2 ) =s (2 ) *const 
g ( 1, 3) =s (3) *const 
g (2, 1) =s (4) *const 
g (2, 2) =s (5) *const 
g (2 , 3) =s ( 6) *const 
g (3 f 1) =s (7) *const 
g (3 f 2) =s (8) *const 
g (3, 3) =s ( 9) *const 
c 

p (1) -s (10) *rmass 
p (2) =s (11) *rmass 
p (3) =s (12) *rmass 
c 

c calculate the tensor h 
c 

do 20 i-1, 3 
do 21 j-1,3 
h(i, j)-0.0 
do 22 13=1,3 

h(i,j)=h(i,j)+g(i,ia)*rj (ia,j) 

22 continue 
21 continue 
2D continue 
c 

return 

end 

c 

subroutine ifdcrk (ido, neq, t , tend, tol, param, y) 
c 

dimension y (33) , yprime (33) , param (50) 

dimension f 1 (33) , f2 (33) , f3 (33) , f4 (33) ,yl (33) ,y2 (33) ,y3 (33) , 
1 y4 (33) 
c 

c external fen 

c 

nstep=1000 
rnstep=nstep 
delt= (tend-t) /rnstep 
c 

do 10 i=l,nstep 
c 

ri=i 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


tl=t+ (ri-1 .0) Melt 
t2=tl+delt/2 .0 
t3=tl+delt/2 .0 
t4=t 1+delt 

do 21 j=l,neq 

yi ( j)=y( j) 

21 continue 

call fen (neq, tl, yl, f 1) 
do 22 j=l,neq 

y2 ( j)«y ( j)+(delt/2 .0) *fl ( j) 

22 continue 

call fen (neq, t2 , y2, f 2) 
do 23 j = l, neq 

y3 ( j ) =y < j ) + (delt /2 . 0 ) *f 2 ( j ) 

23 continue 

call fen (neq, t3, y3 f f 3) 

do 24 j=l , neq 

y4 ( j) =y ( j) +delt*f3 ( j) 

24 continue 

call fen (neq, t 4, y4, f 4) 
do 30 j=l,neq 

y ( j) *y ( j) + (delt/ 6 . 0 ) * <f 1 ( j ) +2 . 0*f2 ( j ) +2 . 0*f 3 ( j ) +f 4 ( j) ) 

30 continue 

10 continue 

return 

end 

subroutine ivderk (ido, neq, tbeg, tend, tol, param, y) 
common/intpar/epsint , ymxref 

dimension y(33) ,yw(33) , ymax (33) , param (50) , err (33) 
dimension fl(33) , f 2 (33) , f 3 (33) , f 4 (33) ,yl<33> ,y2 (33) ,y3 (33) , 
1 y4 (33) 

eps=l . Oe-6 
ymxref =1 . 0e-6 

eps=epsint 

delmin=l . 0e-24 
nstep=1000000 

rnstep=nstep 

delt= (tend-tbeg) /rnstep 

delmin=delt/10 . 0 


c 



fur V 


G-42 



do 11 i=l,neq 

ymax { i ) =amaxl (ymxref , abs (y (i) ) ) 

11 continue 
c 

t=tbeg 

c 

100 tl=t 

if ( t+delt . gt . tend) delt=tend-t 
t2=t+delt /2 . 0 
t3=t+delt 
c 

do 91 i-l,neq 
yl (i) =y (i) 

91 continue 
c 

call rkstep (ido, neq, t 1 , t3, tol, param, yl, y2) 
call rkstep (ido, neq, tl, t2, tol, param, yl, y3) 
call rkstep ( ido, neq, t2 , t 3, tol, param, y3, y4) 
c 

errmax=0 . 0 
do 12 i=l,neq 

err (i) =abs (y4 (i) -y2 (i) ) / (eps*ymax(i) ) 
errmax=amaxl (errmax f err(i) ) 

12 continue 
. c 

if (errmax . le . 0 . 0 ) delt=2 . 0*delt 
if (errmax . gt . 1 . 0 ) go to 101_ 

if {errmax. gt .0.0) delt=0 . 99‘*delt*errmax** {-1 . 0/5 . 0) 
c 

do 21 i=l,neq 
y (i) =y4 (i) 

ymax (i) =amaxl (ymax (i) , abs (y (i) ) ) 

21 continue 
c 

t=t 3 

if (t .ge . tend) go to 99 

if (delt . It . delmin) go to 98 
go to 100 
c 

>101 delt=delt/2.0 

if (delt .It .delmin) go to 98 
go to 100 
c 

98 ido=98 

write { 10 , 110) ido, t, delt 
110 format (lx, i5, 2el5 . 3) 

c 110 format (lx, ’time step less than minimum, ido = T ,i5) 
return 
c 

99 ido=99 

c write ( 10, 120 ) ido, t, delt 

c 120 f ormat ( lx, i5, 2el5 . 3) 
return 
c 

end 

c 

subroutine rkstep (ido, neq, tbeg, tend, tol, param, y, yout) 
c 
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dimension y (33) , pa ram (50) , yout (33) 

dimension f 1 (33) , f 2 (33) , f 3 (33) , f 4 (33) , yl (33) , y2 (33) , y3 (33) , 
1 y4(33) 
c 

delt=tend-tbeg 
t l=tbeg 
t2=t 1+delt /2 . 0 
t3=tl+delt/2 . 0 
t4=tl+delt 
c 

do 21 j=l,neq 

yKj)-y(j) 

21 continue 
c 

call fen (neq, tl, yl, f 1) 
c 

do 22 j=l,neq 

y2 < j) =y ( j ) + (delt/2 . 0) *fl ( j) 

22 continue 
c 

call fen (neq, t2, y2, f2) 
c 

do 23 j=*l, neq 

y3 { j ) =y ( j ) + ( delt/2 . 0 ) *f 2 (j) 

23 continue 
c 

call fen (neq, t3, y3, f 3) 
c 

do 24 j=l,neq 

y4 ( j ) *=y ( j) +delt*f 3 ( j) 

24 continue 
c 

call f cn (neq, t4, y4, f 4) 
c 

do 30 j-^neq 

yout ( j)=y(j) + (delt/6.0)*(fl< j) +2 . 0*f2 ( j) +2 .0*f3 ( j) +f 4 ( j) ) 

30 continue 
c 

return 

end 

c 

subroutine matvl2 (if lag, n, rm, rmi) 
c 

c gauss-jordon elimination 

c 

dimension rm(12, 12) , rmw (12,12), rmi (12,12), chk (12,12), 

1 v ( 12 ) , vw ( 12 ) , it rak ( 12 ) 
c 

c open (1, f ile= T gauss . inp T ) 

c open (2, file- * gauss . out 1 ) 

c 

c write (10 , 101) n 

c 101 format(/*n -*,i5) 
c 

c do 40 i=l,n 

c do 41 j=l,n 

c write (10, 422) rm(i,j) 

c 422 f ormat (el5 . 3) 
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c 41 continue 
c 40 continue 
c 

c do 42 i=l,n 

c read (1, 423) v(i) 

c 423 format (e!5.3) 
c vw ( i ) = v ( i ) 

c 42 continue 
c 

do 20 i=l f n 
c itrak(i)=i 

do 21 j = l,n 
rmw(i, j) =rm ( i, j) 
chk (i, j ) =0 . 0 
rmi (i, j ) =0 . 0 
if (i.eq. j) rmi <i, j) *1 .0 
21 continue 
20 continue 
c 

do 10 i=l,n 
c 

isav=i 

pmag=abs ( rmw ( i , i ) ) 
do 301 irow=i,n 
tmpnum=abs (rmw (irow, i) ) 
if (tmpnum.gt .pmag) isav=irow 

301 continue 

c itrak (i) =isav 

do 302 jcol=l,n 

hold=rmw (i, jcol) 

rmw (i f jcol) =rmw(isav f jcol) 

rmw (isav, jcol) =hold 

hold=rmi (i f jcol) 

rmi (i r jcol) =rmi (isav f jcol) 

rmi (isav, jcol) =hold 

302 continue 
c 

pivot = rmw (i, i) 
c vw (i) =vw (i) /pivot 

'c 

do 11 jj*l,n 

rmi (i, j j) =rmi (i, j j) /pivot 
rmw (if j j)=rmw(i, j j) /pivot 

11 continue 
c 

do 12 k=l , n 
if (k.eq.i) go to 12 
scale=rmw (k, i) 

c vw (k) =vw (k) -scale*vw (i) 

do 13 j-l r n 

rmi (k, j) =rmi (k, j) -scale*rmi (i, j) 
rmw (k, j ) -rmw (k, j ) -scale*rmw (i, j) 

13 continue 

12 continue 
c 

10 continue 
c 

go to 99 
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chksum=0 . 0 
do 50 i=l,n 
do 51 j = 1 , n 
do 52 k=l,n 

chk(i, j ) -chk (i, j)+rm(i,k) *rmi(k, j) 

52 continue 

if (i .eq. j ) chk <i, j ) =chk (i f j) -1 . 0 
chksum^chksum+chk (i, j ) **2 
51 continue 
50 continue 
c 

if (iflag.ne.4) go to 99 
ichk=0 

if (chksum.ge . 1 . 0e-12) ichk=-999 
write (4, 222) ichk, if lag, n, chksum 

222 format (Ix,3i5,el5.3) 
c 

do 30 i=l,n 

do 31 j=l,n 

write (4,223) i,j, rm(i, j ) , rmw (i, j ) , rmi (i, j) , chk (i, j) 

223 format ( lx, 2i5, 4el5 . 3) 

31 continue 

30 continue 
c 

99 dummy=1.0 
c 

c do 44 i-l r n 

c write (2, 223) i,v(i),vw(i) 

c 223 format (lx, i5, 2el5 . 3) 
c 44 continue 
c 

return 

end 

c 

subroutine matv03 (if lag, n, rm, rmi) 
c 

c gauss-jordon elimination 

c 

dimension rm (3, 3) , rmw (3, 3) , rmi (3,3) , chk (3,3) , 

1 v (3) , vw (3) , itrak (3) 
c 

c open (1, file* ’gauss . inp T ) 

c open (2, file* ’gauss . out T ) 

c 

c write ( 10, 101) n 

c 101 format(/ , n * f ,i5) 

c 

c do 40 i=l,n 

c do 41 j=l,n 

c write (10, 422) rm(i,j) 

c 422 f ormat (el5 . 3) 
c 41 continue 
c 40 continue 
c 

c do 42 i=l,n 

c read (1, 423) v (i) 

c 423 format (el5 . 3) 
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c vw(i)=v(i) 

c 42 continue 
c 

do 20 i=l r n 
c itrak(i)=i 

do 21 j-l,n 
rmw (i, j) =rm(i, j) 
chk ( i , j) =0 . 0 
rmi (i r j ) =0 . 0 
if ( i . eq. j ) rmi (i, j ) =1 . 0 
21 continue 
20 continue 
c 

do 10 i=l,n 
c 

isav=i 

pmag=abs (rmw (i, i) ) 

do 301 irow=i,n 

tmpnum=abs (rmw(irow, i) ) 

if ( tmpnum . gt . pmag ) isav=irow 

301 continue 

c itrak (i) =isav 

do 302 jcol=l, n 

hold=rmw (i, jcol) 

rmw (i, jcol) =rmw(isav / jcol) 

rmw (isav, jcol) =hold 

hold=rmi ( i f jcol) 

rmi (i, jcol) =rmi (isav, jcol) 

rmi (isav, jcol) =hold 

302 continue 
c 

pivot=rmw (i, i) 
c vw ( i ) -vw ( i ) /pivot 

c 

do 11 jj=l,n 

rmi (i, j j) =rmi (i, j j) /pivot 
rmw (i, j j)=rmw(i f j j) /pivot 

11 continue 
c 

do 12 k=l , n 
if(k.eq.i) go to 12 
scale=rmw (k, i) 

c vw (k) =vw (k) -scale * vw (i) 

do 13 j=l , n 

rmi (k, j) =rmi (k, j) -scale*rmi <i, j) 
rmw (k, j ) =rmw (k, j) -scale*rmw ( i, j) 

13 continue 

12 continue 
c 

10 continue 
c 

go to 99 

chksum=0 . 0 
do 50 i=l,n 
do 51 j=l,n 
do 52 k=l,n 

chk (i, j) =chk<i, j) +rm(i, k) *rmi (k, j) 


c 
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52 continue 

if (i .eq. j) chk (i, j) =chk (i, j) -1 . 0 
chksum=chksum+chk (i, j) **2 
51 continue 
50 continue 
c 

if ( if lag . ne . 4 ) go to 99 
ichk=0 

if (chksum.ge.l.0e-12) ichk=-999 
write ( 4, 222 ) ichk, if lag, n, chksum 

222 format (Ix,3i5,el5. 3) 
c 

do 30 i=l / n 

do 31 j=l,n 

write (4,223) i,j,rm(i,j), rmw (i, j) , rmi (i, j ) , chk (i, j ) 

223 format (lx,2i5, 4el5.3) 

31 continue 

30 continue 
c 

99 dummy =1 . 0 
c 

c do 44 i-l,n 

c write(2,223) i,v(i),vw(i) 

c 223 format (lx, i5, 2el5 . 3) 
c 44 continue 
c 

return 

end 


u 
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